Differential geometry based multiscale models

Multiscale differential geometry learning of networks with applications to single-cell RNA sequencing data

Single-cell RNA sequencing (scRNA-seq) has emerged as a transformative technology, offering unparalleled insights into the intricate landscape of cellular diversity and gene expression dynamics. scRNA-seq analysis represents a challenging and cutting-edge frontier within the field of biological research. Differential geometry serves as a powerful mathematical tool in various applications of scientific research. In this study, we introduce, for the first time, a multiscale differential geometry (MDG) strategy for addressing the challenges encountered in scRNA-seq data analysis. We assume that intrinsic properties of cells lie on a family of low-dimensional manifolds embedded in the high-dimensional space of scRNA-seq data. Multiscale cell-cell interactive manifolds are constructed to reveal complex relationships in the cell-cell network, where curvature-based features for cells can decipher the intricate structural and biological information. We showcase the utility of our novel approach by demonstrating its effectiveness in classifying cell types. This innovative application of differential geometry in scRNA-seq analysis opens new avenues for understanding the intricacies of biological networks and holds great potential for network analysis in other fields.

Calculation of electrostatic free energy for the nonlinear Poisson-Boltzmann model based on the dimensionless potential

The Poisson-Boltzmann (PB) equation governing the electrostatic potential with a unit is often transformed to a normalized form for a dimensionless potential in numerical studies. To calculate the electrostatic free energy (EFE) of biological interests, a unit conversion has to be conducted, because the existing PB energy functionals are all described in terms of the original potential. To bypass this conversion, this paper proposes energy functionals in terms of the dimensionless potential for the first time in the literature, so that the normalized PB equation can be directly derived by using the Euler-Lagrange variational analysis. Moreover, alternative energy forms have been rigorously derived to avoid approximating the gradient of singular functions in the electrostatic stress term. A systematic study has been carried out to examine the surface integrals involved in alternative energy forms and their dependence on finite domain size and mesh step size, which leads to a recommendation on the EFE forms for efficient computation of protein systems. The calculation of the EFE in the regularization formulation, which is an analytical approach for treating singular charge sources of the PB equation, has also been studied. The proposed energy forms have been validated by considering smooth dielectric settings, such as diffuse interface and super-Gaussian, for which the EFE of the nonlinear PB model is found to be significantly different from that of the linearized PB model. All proposed energy functionals and EFE forms are designed such that the dimensionless potential can be simply plugged in to compute the EFE in the unit of kcal/mol, and they can also be applied in the classical sharp interface PB model.

Finite ion size effects on I-V relations via Poisson-Nernst-Planck systems with two cations: A case study

We consider a quasi-one-dimensional Poisson-Nernst-Planck model with two cations having the same valances and one anion. Bikerman's local hard-sphere potential is included to account for ion size effects. Under some further restrictions on the boundary conditions of the two cations, we obtain approximations of the I-V (current-voltage) relations by treating the ion sizes as small parameters. Critical potentials are identified, which play critical roles in characterizing finite ion size effects on ionic flows. Nonlinear interplays between system parameters, such as boundary concentrations and diffusion coefficients, are analyzed. To provide more intuitive illustrations of our analytical results and better understanding of the dynamics of ionic flows through membrane channels, numerical simulations are performed.

Methodology-Centered Review of Molecular Modeling, Simulation, and Prediction of SARS-CoV-2

Despite tremendous efforts in the past two years, our understanding of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), virus-host interactions, immune response, virulence, transmission, and evolution is still very limited. This limitation calls for further in-depth investigation. Computational studies have become an indispensable component in combating coronavirus disease 2019 (COVID-19) due to their low cost, their efficiency, and the fact that they are free from safety and ethical constraints. Additionally, the mechanism that governs the global evolution and transmission of SARS-CoV-2 cannot be revealed from individual experiments and was discovered by integrating genotyping of massive viral sequences, biophysical modeling of protein-protein interactions, deep mutational data, deep learning, and advanced mathematics. There exists a tsunami of literature on the molecular modeling, simulations, and predictions of SARS-CoV-2 and related developments of drugs, vaccines, antibodies, and diagnostics. To provide readers with a quick update about this literature, we present a comprehensive and systematic methodology-centered review. Aspects such as molecular biophysics, bioinformatics, cheminformatics, machine learning, and mathematics are discussed. This review will be beneficial to researchers who are looking for ways to contribute to SARS-CoV-2 studies and those who are interested in the status of the field.

Geometric algebra generation of molecular surfaces

Geometric algebra is a powerful framework that unifies mathematics and physics. Since its revival in the 1960s, it has attracted great attention and has been exploited in fields like physics, computer science and engineering. This work introduces a geometric algebra method for the molecular surface generation that uses the Clifford-Fourier transform (CFT) which is a generalization of the classical Fourier transform. Notably, the classical Fourier transform and CFT differ in the derivative property in [Formula: see text] for k even. This distinction is due to the non-commutativity of geometric product of pseudoscalars with multivectors and has significant consequences in applications. We use the CFT in [Formula: see text] to benefit from the derivative property in solving partial differential equations (PDEs). The CFT is used to solve the mode decomposition process in PDE transform. Two different initial cases are proposed to make the initial shapes in the present method. The proposed method is applied first to small molecules and proteins. To validate the method, the molecular surfaces generated are compared to surfaces of other definitions. Applications are considered to protein electrostatic surface potentials and solvation free energy. This work opens the door for further applications of geometric algebra and CFT in biological sciences.

EVOLUTIONARY DE RHAM-HODGE METHOD

The de Rham-Hodge theory is a landmark of the 20^th Century's mathematics and has had a great impact on mathematics, physics, computer science, and engineering. This work introduces an evolutionary de Rham-Hodge method to provide a unified paradigm for the multiscale geometric and topological analysis of evolving manifolds constructed from a filtration, which induces a family of evolutionary de Rham complexes. While the present method can be easily applied to close manifolds, the emphasis is given to more challenging compact manifolds with 2-manifold boundaries, which require appropriate analysis and treatment of boundary conditions on differential forms to maintain proper topological properties. Three sets of unique evolutionary Hodge Laplacians are proposed to generate three sets of topology-preserving singular spectra, for which the multiplicities of zero eigenvalues correspond to exactly the persistent Betti numbers of dimensions 0, 1 and 2. Additionally, three sets of non-zero eigenvalues further reveal both topological persistence and geometric progression during the manifold evolution. Extensive numerical experiments are carried out via the discrete exterior calculus to demonstrate the potential of the proposed paradigm for data representation and shape analysis of both point cloud data and density maps. To demonstrate the utility of the proposed method, the application is considered to the protein B-factor predictions of a few challenging cases for which existing biophysical models break down.

Competition between Cations via Classical Poisson-Nernst-Planck Models with Nonzero but Small Permanent Charges

We study a one-dimensional Poisson–Nernst–Planck system for ionic flow through a membrane channel. Nonzero but small permanent charge, the major structural quantity of an ion channel, is included in the model. Two cations with the same valences and one anion are included in the model, which provides more rich and complicated correlations/interactions between ions. The cross-section area of the channel is included in the system, and it provides certain information of the geometry of the three-dimensional channel, which is critical for our analysis. Geometric singular perturbation analysis is employed to establish the existence and local uniqueness of solutions to the system for small permanent charges. Treating the permanent charge as a small parameter, through regular perturbation analysis, we are able to derive approximations of the individual fluxes explicitly, and this allows us to study the competition between two cations, which is related to the selectivity phenomena of ion channels. Numerical simulations are performed to provide a more intuitive illustration of our analytical results, and they are consistent.

A Computational Model of Protein Induced Membrane Morphology with Geodesic Curvature Driven Protein-Membrane Interface

Continuum or hybrid modeling of bilayer membrane morphological dynamics induced by embedded proteins necessitates the identification of protein-membrane interfaces and coupling of deformations of two surfaces. In this article we developed (i) a minimal total geodesic curvature model to describe these interfaces, and (ii) a numerical one-one mapping between two surface through a conformal mapping of each surface to the common middle annulus. Our work provides the first computational tractable approach for determining the interfaces between bilayer and embedded proteins. The one-one mapping allows a convenient coupling of the morphology of two surfaces. We integrated these two new developments into the energetic model of protein-membrane interactions, and developed the full set of numerical methods for the coupled system. Numerical examples are presented to demonstrate (1) the efficiency and robustness of our methods in locating the curves with minimal total geodesic curvature on highly complicated protein surfaces, (2) the usefulness of these interfaces as interior boundaries for membrane deformation, and (3) the rich morphology of bilayer surfaces for different protein-membrane interfaces.

A review of mathematical representations of biomolecular data

Recently, machine learning (ML) has established itself in various worldwide benchmarking competitions in computational biology, including Critical Assessment of Structure Prediction (CASP) and Drug Design Data Resource (D3R) Grand Challenges. However, the intricate structural complexity and high ML dimensionality of biomolecular datasets obstruct the efficient application of ML algorithms in the field. In addition to data and algorithm, an efficient ML machinery for biomolecular predictions must include structural representation as an indispensable component. Mathematical representations that simplify the biomolecular structural complexity and reduce ML dimensionality have emerged as a prime winner in D3R Grand Challenges. This review is devoted to the recent advances in developing low-dimensional and scalable mathematical representations of biomolecules in our laboratory. We discuss three classes of mathematical approaches, including algebraic topology, differential geometry, and graph theory. We elucidate how the physical and biological challenges have guided the evolution and development of these mathematical apparatuses for massive and diverse biomolecular data. We focus the performance analysis on protein-ligand binding predictions in this review although these methods have had tremendous success in many other applications, such as protein classification, virtual screening, and the predictions of solubility, solvation free energies, toxicity, partition coefficients, protein folding stability changes upon mutation, etc.

MathDL: mathematical deep learning for D3R Grand Challenge 4

We present the performances of our mathematical deep learning (MathDL) models for D3R Grand Challenge 4 (GC4). This challenge involves pose prediction, affinity ranking, and free energy estimation for beta secretase 1 (BACE) as well as affinity ranking and free energy estimation for Cathepsin S (CatS). We have developed advanced mathematics, namely differential geometry, algebraic graph, and/or algebraic topology, to accurately and efficiently encode high dimensional physical/chemical interactions into scalable low-dimensional rotational and translational invariant representations. These representations are integrated with deep learning models, such as generative adversarial networks (GAN) and convolutional neural networks (CNN) for pose prediction and energy evaluation, respectively. Overall, our MathDL models achieved the top place in pose prediction for BACE ligands in Stage 1a. Moreover, our submissions obtained the highest Spearman correlation coefficient on the affinity ranking of 460 CatS compounds, and the smallest centered root mean square error on the free energy set of 39 CatS molecules. It is worthy to mention that our method on docking pose predictions has significantly improved from our previous ones.

Nonlinear simulation of an elastic tumor-host interface

We develop a computational method for simulating the nonlinear dynamics of an elastic tumor-host interface. This work is motivated by the recent linear stability analysis of a two-phase tumor model with an elastic membrane interface in 2D [47]. Unlike the classic tumor model with surface tension, the elastic interface condition is numerically challenging due to the 4th order derivative from the Helfrich bending energy. Here we are interested in exploring the nonlinear interface dynamics in a sharp interface framework. We consider a curvature dependent bending rigidity (curvature weakening [22]) to investigate metastasis patterns such as chains or fingers that invade the host environment. We solve the nutrient field and the Stokes flow field using a spectrally accurate boundary integral method, and update the interface using a nonstiff semi-implicit approach. Numerical results suggest curvature weakening promotes the development of branching patterns instead of encapsulated morphologies in a long period of time. For non-weakened bending rigidity, we are able to find self-similar shrinking morphologies based on marginally stable value of the apoptosis rate.

DG-GL: Differential geometry-based geometric learning of molecular datasets

MOTIVATION

Despite its great success in various physical modeling, differential geometry (DG) has rarely been devised as a versatile tool for analyzing large, diverse, and complex molecular and biomolecular datasets because of the limited understanding of its potential power in dimensionality reduction and its ability to encode essential chemical and biological information in differentiable manifolds.

RESULTS

We put forward a differential geometry-based geometric learning (DG-GL) hypothesis that the intrinsic physics of three-dimensional (3D) molecular structures lies on a family of low-dimensional manifolds embedded in a high-dimensional data space. We encode crucial chemical, physical, and biological information into 2D element interactive manifolds, extracted from a high-dimensional structural data space via a multiscale discrete-to-continuum mapping using differentiable density estimators. Differential geometry apparatuses are utilized to construct element interactive curvatures in analytical forms for certain analytically differentiable density estimators. These low-dimensional differential geometry representations are paired with a robust machine learning algorithm to showcase their descriptive and predictive powers for large, diverse, and complex molecular and biomolecular datasets. Extensive numerical experiments are carried out to demonstrate that the proposed DG-GL strategy outperforms other advanced methods in the predictions of drug discovery-related protein-ligand binding affinity, drug toxicity, and molecular solvation free energy.

AVAILABILITY AND IMPLEMENTATION

http://weilab.math.msu.edu/DG-GL/ Contact: wei@math.msu.edu.

Generative network complex (GNC) for drug discovery

It remains a challenging task to generate a vast variety of novel compounds with desirable pharmacological properties. In this work, a generative network complex (GNC) is proposed as a new platform for designing novel compounds, predicting their physical and chemical properties, and selecting potential drug candidates that fulfill various druggable criteria such as binding affinity, solubility, partition coefficient, etc. We combine a SMILES string generator, which consists of an encoder, a drug-property controlled or regulated latent space, and a decoder, with verification deep neural networks, a target-specific three-dimensional (3D) pose generator, and mathematical deep learning networks to generate new compounds, predict their drug properties, construct 3D poses associated with target proteins, and reevaluate druggability, respectively. New compounds were generated in the latent space by either randomized output, controlled output, or optimized output. In our demonstration, 2.08 million and 2.8 million novel compounds are generated respectively for Cathepsin S and BACE targets. These new compounds are very different from the seeds and cover a larger chemical space. For potentially active compounds, their 3D poses are generated using a state-of-the-art method. The resulting 3D complexes are further evaluated for druggability by a championing deep learning algorithm based on algebraic topology, differential geometry, and algebraic graph theories. Performed on supercomputers, the whole process took less than one week. Therefore, our GNC is an efficient new paradigm for discovering new drug candidates.

Divide-and-conquer strategy for large-scale Eulerian solvent excluded surface

MOTIVATION

Surface generation and visualization are some of the most important tasks in biomolecular modeling and computation. Eulerian solvent excluded surface (ESES) software provides analytical solvent excluded surface (SES) in the Cartesian grid, which is necessary for simulating many biomolecular electrostatic and ion channel models. However, large biomolecules and/or fine grid resolutions give rise to excessively large memory requirements in ESES construction. We introduce an out-of-core and parallel algorithm to improve the ESES software.

RESULTS

The present approach drastically improves the spatial and temporal efficiency of ESES. The memory footprint and time complexity are analyzed and empirically verified through extensive tests with a large collection of biomolecule examples. Our results show that our algorithm can successfully reduce memory footprint through a straightforward divide-and-conquer strategy to perform the calculation of arbitrarily large proteins on a typical commodity personal computer. On multi-core computers or clusters, our algorithm can reduce the execution time by parallelizing most of the calculation as disjoint subproblems. Various comparisons with the state-of-the-art Cartesian grid based SES calculation were done to validate the present method and show the improved efficiency. This approach makes ESES a robust software for the construction of analytical solvent excluded surfaces.

AVAILABILITY AND IMPLEMENTATION

http://weilab.math.msu.edu/ESES.

Fractional Poisson-Nernst-Planck Model for Ion Channels I: Basic Formulations and Algorithms

In this work, we propose a fractional Poisson-Nernst-Planck model to describe ion permeation in gated ion channels. Due to the intrinsic conformational changes, crowdedness in narrow channel pores, binding and trapping introduced by functioning units of channel proteins, ionic transport in the channel exhibits a power-law-like anomalous diffusion dynamics. We start from continuous-time random walk model for a single ion and use a long-tailed density distribution function for the particle jump waiting time, to derive the fractional Fokker-Planck equation. Then, it is generalized to the macroscopic fractional Poisson-Nernst-Planck model for ionic concentrations. Necessary computational algorithms are designed to implement numerical simulations for the proposed model, and the dynamics of gating current is investigated. Numerical simulations show that the fractional PNP model provides a more qualitatively reasonable match to the profile of gating currents from experimental observations. Meanwhile, the proposed model motivates new challenges in terms of mathematical modeling and computations.

Generalized flexibility-rigidity index

Flexibility-rigidity index (FRI) has been developed as a robust, accurate, and efficient method for macromolecular thermal fluctuation analysis and B-factor prediction. The performance of FRI depends on its formulations of rigidity index and flexibility index. In this work, we introduce alternative rigidity and flexibility formulations. The structure of the classic Gaussian surface is utilized to construct a new type of rigidity index, which leads to a new class of rigidity densities with the classic Gaussian surface as a special case. Additionally, we introduce a new type of flexibility index based on the domain indicator property of normalized rigidity density. These generalized FRI (gFRI) methods have been extensively validated by the B-factor predictions of 364 proteins. Significantly outperforming the classic Gaussian network model, gFRI is a new generation of methodologies for accurate, robust, and efficient analysis of protein flexibility and fluctuation. Finally, gFRI based molecular surface generation and flexibility visualization are demonstrated.

ESES: Software for Eulerian solvent excluded surface

Solvent excluded surface (SES) is one of the most popular surface definitions in biophysics and molecular biology. In addition to its usage in biomolecular visualization, it has been widely used in implicit solvent models, in which SES is usually immersed in a Cartesian mesh. Therefore, it is important to construct SESs in the Eulerian representation for biophysical modeling and computation. This work describes a software package called Eulerian solvent excluded surface (ESES) for the generation of accurate SESs in Cartesian grids. ESES offers the description of the solvent and solute domains by specifying all the intersection points between the SES and the Cartesian grid lines. Additionally, the interface normal at each intersection point is evaluated. Furthermore, for a given biomolecule, the ESES software not only provides the whole surface area, but also partitions the surface area according to atomic types. Homology theory is utilized to detect topological features, such as loops and cavities, on the complex formed by the SES. The sizes of loops and cavities are measured based on persistent homology with an evolutionary partial differential equation-based filtration. ESES is extensively validated by surface visualization, electrostatic solvation free energy computation, surface area and volume calculations, and loop and cavity detection and their size estimation. We used the Amber PBSA test set in our electrostatic solvation energy, area, and volume validations. Our results are either calibrated by analytical values or compared with those from the MSMS software. © 2017 Wiley Periodicals, Inc.

Fast Simulation of Lipid Vesicle Deformation Using Spherical Harmonic Approximation

Lipid vesicles appear ubiquitously in biological systems. Understanding how the mechanical and intermolecular interactions deform vesicle membranes is a fundamental question in biophysics. In this article we develop a fast algorithm to compute the surface configurations of lipid vesicles by introducing surface harmonic functions to approximate the membrane surface. This parameterization allows an analytical computation of the membrane curvature energy and its gradient for the efficient minimization of the curvature energy using a nonlinear conjugate gradient method. Our approach drastically reduces the degrees of freedom for approximating the membrane surfaces compared to the previously developed finite element and finite difference methods. Vesicle deformations with a reduced volume larger than 0.65 can be well approximated by using as small as 49 surface harmonic functions. The method thus has a great potential to reduce the computational expense of tracking multiple vesicles which deform for their interaction with external fields.

The impact of surface area, volume, curvature, and Lennard-Jones potential to solvation modeling

This article explores the impact of surface area, volume, curvature, and Lennard-Jones (LJ) potential on solvation free energy predictions. Rigidity surfaces are utilized to generate robust analytical expressions for maximum, minimum, mean, and Gaussian curvatures of solvent-solute interfaces, and define a generalized Poisson-Boltzmann (GPB) equation with a smooth dielectric profile. Extensive correlation analysis is performed to examine the linear dependence of surface area, surface enclosed volume, maximum curvature, minimum curvature, mean curvature, and Gaussian curvature for solvation modeling. It is found that surface area and surfaces enclosed volumes are highly correlated to each other's, and poorly correlated to various curvatures for six test sets of molecules. Different curvatures are weakly correlated to each other for six test sets of molecules, but are strongly correlated to each other within each test set of molecules. Based on correlation analysis, we construct twenty six nontrivial nonpolar solvation models. Our numerical results reveal that the LJ potential plays a vital role in nonpolar solvation modeling, especially for molecules involving strong van der Waals interactions. It is found that curvatures are at least as important as surface area or surface enclosed volume in nonpolar solvation modeling. In conjugation with the GPB model, various curvature-based nonpolar solvation models are shown to offer some of the best solvation free energy predictions for a wide range of test sets. For example, root mean square errors from a model constituting surface area, volume, mean curvature, and LJ potential are less than 0.42 kcal/mol for all test sets. © 2016 Wiley Periodicals, Inc.

Multiscale method for modeling binding phenomena involving large objects: application to kinesin motor domains motion along microtubules

Many biological phenomena involve the binding of proteins to a large object. Because the electrostatic forces that guide binding act over large distances, truncating the size of the system to facilitate computational modeling frequently yields inaccurate results. Our multiscale approach implements a computational focusing method that permits computation of large systems without truncating the electrostatic potential and achieves the high resolution required for modeling macromolecular interactions, all while keeping the computational time reasonable. We tested our approach on the motility of various kinesin motor domains. We found that electrostatics help guide kinesins as they walk: N-kinesins towards the plus-end, and C-kinesins towards the minus-end of microtubules. Our methodology enables computation in similar, large systems including protein binding to DNA, viruses, and membranes.

Multiscale models of compact bone

This work is concerned about multiscale models of compact bone. We focus on the lacuna–canalicular system. The interstitial fluid and the ions in it are regarded as solvent and others are treated as solute. The system has the characteristic of solvation process as well as non-equilibrium dynamics. The differential geometry theory of surfaces is adopted. We use this theory to separate the macroscopic domain of solvent from the microscopic domain of solute. We also use it to couple continuum and discrete descriptions. The energy functionals are constructed and then the variational principle is applied to the energy functionals so as to derive desirable governing equations. We consider both long-range polar interactions and short-range nonpolar interactions. The solution of governing equations leads to the minimization of the total energy.

Object-oriented Persistent Homology

Persistent homology provides a new approach for the topological simplification of big data via measuring the life time of intrinsic topological features in a filtration process and has found its success in scientific and engineering applications. However, such a success is essentially limited to qualitative data classification and analysis. Indeed, persistent homology has rarely been employed for quantitative modeling and prediction. Additionally, the present persistent homology is a passive tool, rather than a proactive technique, for classification and analysis. In this work, we outline a general protocol to construct object-oriented persistent homology methods. By means of differential geometry theory of surfaces, we construct an objective functional, namely, a surface free energy defined on the data of interest. The minimization of the objective functional leads to a Laplace-Beltrami operator which generates a multiscale representation of the initial data and offers an objective oriented filtration process. The resulting differential geometry based object-oriented persistent homology is able to preserve desirable geometric features in the evolutionary filtration and enhances the corresponding topological persistence. The cubical complex based homology algorithm is employed in the present work to be compatible with the Cartesian representation of the Laplace-Beltrami flow. The proposed Laplace-Beltrami flow based persistent homology method is extensively validated. The consistence between Laplace-Beltrami flow based filtration and Euclidean distance based filtration is confirmed on the Vietoris-Rips complex for a large amount of numerical tests. The convergence and reliability of the present Laplace-Beltrami flow based cubical complex filtration approach are analyzed over various spatial and temporal mesh sizes. The Laplace-Beltrami flow based persistent homology approach is utilized to study the intrinsic topology of proteins and fullerene molecules. Based on a quantitative model which correlates the topological persistence of fullerene central cavity with the total curvature energy of the fullerene structure, the proposed method is used for the prediction of fullerene isomer stability. The efficiency and robustness of the present method are verified by more than 500 fullerene molecules. It is shown that the proposed persistent homology based quantitative model offers good predictions of total curvature energies for ten types of fullerene isomers. The present work offers the first example to design object-oriented persistent homology to enhance or preserve desirable features in the original data during the filtration process and then automatically detect or extract the corresponding topological traits from the data.

Parameter optimization in differential geometry based solvation models

Differential geometry (DG) based solvation models are a new class of variational implicit solvent approaches that are able to avoid unphysical solvent-solute boundary definitions and associated geometric singularities, and dynamically couple polar and non-polar interactions in a self-consistent framework. Our earlier study indicates that DG based non-polar solvation model outperforms other methods in non-polar solvation energy predictions. However, the DG based full solvation model has not shown its superiority in solvation analysis, due to its difficulty in parametrization, which must ensure the stability of the solution of strongly coupled nonlinear Laplace-Beltrami and Poisson-Boltzmann equations. In this work, we introduce new parameter learning algorithms based on perturbation and convex optimization theories to stabilize the numerical solution and thus achieve an optimal parametrization of the DG based solvation models. An interesting feature of the present DG based solvation model is that it provides accurate solvation free energy predictions for both polar and non-polar molecules in a unified formulation. Extensive numerical experiment demonstrates that the present DG based solvation model delivers some of the most accurate predictions of the solvation free energies for a large number of molecules.

Matched Interface and Boundary Method for Elasticity Interface Problems

Elasticity theory is an important component of continuum mechanics and has had widely spread applications in science and engineering. Material interfaces are ubiquity in nature and man-made devices, and often give rise to discontinuous coefficients in the governing elasticity equations. In this work, the matched interface and boundary (MIB) method is developed to address elasticity interface problems. Linear elasticity theory for both isotropic homogeneous and inhomogeneous media is employed. In our approach, Lamé's parameters can have jumps across the interface and are allowed to be position dependent in modeling isotropic inhomogeneous material. Both strong discontinuity, i.e., discontinuous solution, and weak discontinuity, namely, discontinuous derivatives of the solution, are considered in the present study. In the proposed method, fictitious values are utilized so that the standard central finite different schemes can be employed regardless of the interface. Interface jump conditions are enforced on the interface, which in turn, accurately determines fictitious values. We design new MIB schemes to account for complex interface geometries. In particular, the cross derivatives in the elasticity equations are difficult to handle for complex interface geometries. We propose secondary fictitious values and construct geometry based interpolation schemes to overcome this difficulty. Numerous analytical examples are used to validate the accuracy, convergence and robustness of the present MIB method for elasticity interface problems with both small and large curvatures, strong and weak discontinuities, and constant and variable coefficients. Numerical tests indicate second order accuracy in both L∞ and L2 norms.

Persistent topology for cryo-EM data analysis

In this work, we introduce persistent homology for the analysis of cryo-electron microscopy (cryo-EM) density maps. We identify the topological fingerprint or topological signature of noise, which is widespread in cryo-EM data. For low signal-to-noise ratio (SNR) volumetric data, intrinsic topological features of biomolecular structures are indistinguishable from noise. To remove noise, we employ geometric flows that are found to preserve the intrinsic topological fingerprints of cryo-EM structures and diminish the topological signature of noise. In particular, persistent homology enables us to visualize the gradual separation of the topological fingerprints of cryo-EM structures from those of noise during the denoising process, which gives rise to a practical procedure for prescribing a noise threshold to extract cryo-EM structure information from noise contaminated data after certain iterations of the geometric flow equation. To further demonstrate the utility of persistent homology for cryo-EM data analysis, we consider a microtubule intermediate structure Electron Microscopy Data (EMD 1129). Three helix models, an alpha-tubulin monomer model, an alpha-tubulin and beta-tubulin model, and an alpha-tubulin and beta-tubulin dimer model, are constructed to fit the cryo-EM data. The least square fitting leads to similarly high correlation coefficients, which indicates that structure determination via optimization is an ill-posed inverse problem. However, these models have dramatically different topological fingerprints. Especially, linkages or connectivities that discriminate one model from another, play little role in the traditional density fitting or optimization but are very sensitive and crucial to topological fingerprints. The intrinsic topological features of the microtubule data are identified after topological denoising. By a comparison of the topological fingerprints of the original data and those of three models, we found that the third model is topologically favored. The present work offers persistent homology based new strategies for topological denoising and for resolving ill-posed inverse problems.

Second order Method for Solving 3D Elasticity Equations with Complex Interfaces

Elastic materials are ubiquitous in nature and indispensable components in man-made devices and equipments. When a device or equipment involves composite or multiple elastic materials, elasticity interface problems come into play. The solution of three dimensional (3D) elasticity interface problems is significantly more difficult than that of elliptic counterparts due to the coupled vector components and cross derivatives in the governing elasticity equation. This work introduces the matched interface and boundary (MIB) method for solving 3D elasticity interface problems. The proposed MIB elasticity interface scheme utilizes fictitious values on irregular grid points near the material interface to replace function values in the discretization so that the elasticity equation can be discretized using the standard finite difference schemes as if there were no material interface. The interface jump conditions are rigorously enforced on the intersecting points between the interface and the mesh lines. Such an enforcement determines the fictitious values. A number of new techniques has been developed to construct efficient MIB elasticity interface schemes for dealing with cross derivative in coupled governing equations. The proposed method is extensively validated over both weak and strong discontinuity of the solution, both piecewise constant and position-dependent material parameters, both smooth and nonsmooth interface geometries, and both small and large contrasts in the Poisson's ratio and shear modulus across the interface. Numerical experiments indicate that the present MIB method is of second order convergence in both L∞ and L2 error norms for handling arbitrarily complex interfaces, including biomolecular surfaces. To our best knowledge, this is the first elasticity interface method that is able to deliver the second convergence for the molecular surfaces of proteins..

Multidimensional persistence in biomolecular data

Persistent homology has emerged as a popular technique for the topological simplification of big data, including biomolecular data. Multidimensional persistence bears considerable promise to bridge the gap between geometry and topology. However, its practical and robust construction has been a challenge. We introduce two families of multidimensional persistence, namely pseudomultidimensional persistence and multiscale multidimensional persistence. The former is generated via the repeated applications of persistent homology filtration to high-dimensional data, such as results from molecular dynamics or partial differential equations. The latter is constructed via isotropic and anisotropic scales that create new simiplicial complexes and associated topological spaces. The utility, robustness, and efficiency of the proposed topological methods are demonstrated via protein folding, protein flexibility analysis, the topological denoising of cryoelectron microscopy data, and the scale dependence of nanoparticles. Topological transition between partial folded and unfolded proteins has been observed in multidimensional persistence. The separation between noise topological signatures and molecular topological fingerprints is achieved by the Laplace-Beltrami flow. The multiscale multidimensional persistent homology reveals relative local features in Betti-0 invariants and the relatively global characteristics of Betti-1 and Betti-2 invariants.

Persistent homology for the quantitative prediction of fullerene stability

Persistent homology is a relatively new tool often used for qualitative analysis of intrinsic topological features in images and data originated from scientific and engineering applications. In this article, we report novel quantitative predictions of the energy and stability of fullerene molecules, the very first attempt in using persistent homology in this context. The ground-state structures of a series of small fullerene molecules are first investigated with the standard Vietoris-Rips complex. We decipher all the barcodes, including both short-lived local bars and long-lived global bars arising from topological invariants, and associate them with fullerene structural details. Using accumulated bar lengths, we build quantitative models to correlate local and global Betti-2 bars, respectively with the heat of formation and total curvature energies of fullerenes. It is found that the heat of formation energy is related to the local hexagonal cavities of small fullerenes, while the total curvature energies of fullerene isomers are associated with their sphericities, which are measured by the lengths of their long-lived Betti-2 bars. Excellent correlation coefficients (>0.94) between persistent homology predictions and those of quantum or curvature analysis have been observed. A correlation matrix based filtration is introduced to further verify our findings.

Persistent homology analysis of protein structure, flexibility, and folding

Proteins are the most important biomolecules for living organisms. The understanding of protein structure, function, dynamics, and transport is one of the most challenging tasks in biological science. In the present work, persistent homology is, for the first time, introduced for extracting molecular topological fingerprints (MTFs) based on the persistence of molecular topological invariants. MTFs are utilized for protein characterization, identification, and classification. The method of slicing is proposed to track the geometric origin of protein topological invariants. Both all-atom and coarse-grained representations of MTFs are constructed. A new cutoff-like filtration is proposed to shed light on the optimal cutoff distance in elastic network models. On the basis of the correlation between protein compactness, rigidity, and connectivity, we propose an accumulated bar length generated from persistent topological invariants for the quantitative modeling of protein flexibility. To this end, a correlation matrix-based filtration is developed. This approach gives rise to an accurate prediction of the optimal characteristic distance used in protein B-factor analysis. Finally, MTFs are employed to characterize protein topological evolution during protein folding and quantitatively predict the protein folding stability. An excellent consistence between our persistent homology prediction and molecular dynamics simulation is found. This work reveals the topology-function relationship of proteins.

Kinetic model for the large deformation of cylindrical gels

Here we investigate the kinetic model for the large deformation theory of hydrogel under the outside stimulations. We present the large deformation dynamical model in the following two points. (1) The phase transition caused by the deformation gradient is concerned, which makes the model more integral. (2) Based on the steady-state model, the time-dependent large deformation model is proposed and the time developing process is investigated. Considering the force of the large deformation, we introduce the heat equation to express the transformation of the chemical potential. The extended model can be used to describe the development of the large deformation. We present numerical examples of the cylindrical hydrogel for one- and two-dimensional cases under pressure and stretch. Besides, some key parameters are studied to test the model.

A fast alternating direction implicit algorithm for geometric flow equations in biomolecular surface generation

In this paper, a new alternating direction implicit (ADI) method is introduced to solve potential driven geometric flow PDEs for biomolecular surface generation. For such PDEs, an extra factor is usually added to stabilize the explicit time integration. However, two existing implicit ADI schemes are also based on the scaled form, which involves nonlinear cross derivative terms that have to be evaluated explicitly. This affects the stability and accuracy of these ADI schemes. To overcome these difficulties, we propose a new ADI algorithm based on the unscaled form so that cross derivatives are not involved. Central finite differences are employed to discretize the nonhomogenous diffusion process of the geometric flow. The proposed ADI algorithm is validated through benchmark examples with analytical solutions, reference solutions, or literature results. Moreover, quantitative indicators of a biomolecular surface, including surface area, surface-enclosed volume, and solvation free energy, are analyzed for various proteins. The proposed ADI method is found to be unconditionally stable and more accurate than the existing ADI schemes in all tests. This enables the use of a large time increment in the steady state simulation so that the proposed ADI algorithm is very efficient for biomolecular surface generation.

Multiscale geometric modeling of macromolecules I: Cartesian representation

This paper focuses on the geometric modeling and computational algorithm development of biomolecular structures from two data sources: Protein Data Bank (PDB) and Electron Microscopy Data Bank (EMDB) in the Eulerian (or Cartesian) representation. Molecular surface (MS) contains non-smooth geometric singularities, such as cusps, tips and self-intersecting facets, which often lead to computational instabilities in molecular simulations, and violate the physical principle of surface free energy minimization. Variational multiscale surface definitions are proposed based on geometric flows and solvation analysis of biomolecular systems. Our approach leads to geometric and potential driven Laplace-Beltrami flows for biomolecular surface evolution and formation. The resulting surfaces are free of geometric singularities and minimize the total free energy of the biomolecular system. High order partial differential equation (PDE)-based nonlinear filters are employed for EMDB data processing. We show the efficacy of this approach in feature-preserving noise reduction. After the construction of protein multiresolution surfaces, we explore the analysis and characterization of surface morphology by using a variety of curvature definitions. Apart from the classical Gaussian curvature and mean curvature, maximum curvature, minimum curvature, shape index, and curvedness are also applied to macromolecular surface analysis for the first time. Our curvature analysis is uniquely coupled to the analysis of electrostatic surface potential, which is a by-product of our variational multiscale solvation models. As an expository investigation, we particularly emphasize the numerical algorithms and computational protocols for practical applications of the above multiscale geometric models. Such information may otherwise be scattered over the vast literature on this topic. Based on the curvature and electrostatic analysis from our multiresolution surfaces, we introduce a new concept, the polarized curvature, for the prediction of protein binding sites.

Modeling and computation of heterogeneous implicit solvent and its applications for biomolecules

Description of inhomogeneous dielectric properties of a solvent in the vicinity of ions has been attracting research interests in mathematical modeling for many years. From many experimental results, it has been concluded that the dielectric response of a solvent linearly depends on the ionic strength within a certain range. Based on this assumption, a new implicit solvent model is proposed in the form of total free energy functional and a quasi-linear Poisson-Boltzmann equation (QPBE) is derived. Classical Newton’s iteration can be used to solve the QPBE numerically but the corresponding Jacobian matrix is complicated due to the quasi-linear term. In the current work, a systematic formulation of the Jacobian matrix is derived. As an alternative option, an algorithm mixing the Newton’s iteration and the fixed point method is proposed to avoid the complicated Jacobian matrix, and it is a more general algorithm for equation with discontinuous coefficients. Computational efficiency and accuracy for these two methods are investigated based on a set of equation parameters. At last, the QPBE with singular charge source and piece-wisely defined dielectric functions has been applied to analyze electrostatics of macro biomolecules in a complicated solvent. A set of computational algorithms such as interface method, singular charge removal technique and the Newtonfixed- point iteration are employed to solve the QPBE. Biological applications of the proposed model and algorithms are provided, including calculation of electrostatic solvation free energy of proteins, investigation of physical properties of channel pore of an ion channel, and electrostatics analysis for the segment of a DNA strand.

Multiscale Multiphysics and Multidomain Models I: Basic Theory

This work extends our earlier two-domain formulation of a differential geometry based multiscale paradigm into a multidomain theory, which endows us the ability to simultaneously accommodate multiphysical descriptions of aqueous chemical, physical and biological systems, such as fuel cells, solar cells, nanofluidics, ion channels, viruses, RNA polymerases, molecular motors and large macromolecular complexes. The essential idea is to make use of the differential geometry theory of surfaces as a natural means to geometrically separate the macroscopic domain of solvent from the microscopic domain of solute, and dynamically couple continuum and discrete descriptions. Our main strategy is to construct energy functionals to put on an equal footing of multiphysics, including polar (i.e., electrostatic) solvation, nonpolar solvation, chemical potential, quantum mechanics, fluid mechanics, molecular mechanics, coarse grained dynamics and elastic dynamics. The variational principle is applied to the energy functionals to derive desirable governing equations, such as multidomain Laplace-Beltrami (LB) equations for macromolecular morphologies, multidomain Poisson-Boltzmann (PB) equation or Poisson equation for electrostatic potential, generalized Nernst-Planck (NP) equations for the dynamics of charged solvent species, generalized Navier-Stokes (NS) equation for fluid dynamics, generalized Newton's equations for molecular dynamics (MD) or coarse-grained dynamics and equation of motion for elastic dynamics. Unlike the classical PB equation, our PB equation is an integral-differential equation due to solvent-solute interactions. To illustrate the proposed formalism, we have explicitly constructed three models, a multidomain solvation model, a multidomain charge transport model and a multidomain chemo-electro-fluid-MD-elastic model. Each solute domain is equipped with distinct surface tension, pressure, dielectric function, and charge density distribution. In addition to long-range Coulombic interactions, various non-electrostatic solvent-solute interactions are considered in the present modeling. We demonstrate the consistency between the non-equilibrium charge transport model and the equilibrium solvation model by showing the systematical reduction of the former to the latter at equilibrium. This paper also offers a brief review of the field.

OPTIMIZATION BIAS IN ENERGY-BASED STRUCTURE PREDICTION

Physics-based computational approaches to predicting the structure of macromolecules such as proteins are gaining increased use, but there are remaining challenges. In the current work, it is demonstrated that in energy-based prediction methods, the degree of optimization of the sampled structures can influence the prediction results. In particular, discrepancies in the degree of local sampling can bias the predictions in favor of the oversampled structures by shifting the local probability distributions of the minimum sampled energies. In simple systems, it is shown that the magnitude of the errors can be calculated from the energy surface, and for certain model systems, derived analytically. Further, it is shown that for energy wells whose forms differ only by a randomly assigned energy shift, the optimal accuracy of prediction is achieved when the sampling around each structure is equal. Energy correction terms can be used in cases of unequal sampling to reproduce the total probabilities that would occur under equal sampling, but optimal corrections only partially restore the prediction accuracy lost to unequal sampling. For multiwell systems, the determination of the correction terms is a multibody problem; it is shown that the involved cross-correlation multiple integrals can be reduced to simpler integrals. The possible implications of the current analysis for macromolecular structure prediction are discussed.

Multiscale multiphysics and multidomain models--flexibility and rigidity

The emerging complexity of large macromolecules has led to challenges in their full scale theoretical description and computer simulation. Multiscale multiphysics and multidomain models have been introduced to reduce the number of degrees of freedom while maintaining modeling accuracy and achieving computational efficiency. A total energy functional is constructed to put energies for polar and nonpolar solvation, chemical potential, fluid flow, molecular mechanics, and elastic dynamics on an equal footing. The variational principle is utilized to derive coupled governing equations for the above mentioned multiphysical descriptions. Among these governing equations is the Poisson-Boltzmann equation which describes continuum electrostatics with atomic charges. The present work introduces the theory of continuum elasticity with atomic rigidity (CEWAR). The essence of CEWAR is to formulate the shear modulus as a continuous function of atomic rigidity. As a result, the dynamics complexity of a macromolecular system is separated from its static complexity so that the more time-consuming dynamics is handled with continuum elasticity theory, while the less time-consuming static analysis is pursued with atomic approaches. We propose a simple method, flexibility-rigidity index (FRI), to analyze macromolecular flexibility and rigidity in atomic detail. The construction of FRI relies on the fundamental assumption that protein functions, such as flexibility, rigidity, and energy, are entirely determined by the structure of the protein and its environment, although the structure is in turn determined by all the interactions. As such, the FRI measures the topological connectivity of protein atoms or residues and characterizes the geometric compactness of the protein structure. As a consequence, the FRI does not resort to the interaction Hamiltonian and bypasses matrix diagonalization, which underpins most other flexibility analysis methods. FRI's computational complexity is of O(N(2)) at most, where N is the number of atoms or residues, in contrast to O(N(3)) for Hamiltonian based methods. We demonstrate that the proposed FRI gives rise to accurate prediction of protein B-Factor for a set of 263 proteins. We show that a parameter free FRI is able to achieve about 95% accuracy of the parameter optimized FRI. An interpolation algorithm is developed to construct continuous atomic flexibility functions for visualization and use with CEWAR.

Multiscale geometric modeling of macromolecules II: Lagrangian representation

Geometric modeling of biomolecules plays an essential role in the conceptualization of biolmolecular structure, function, dynamics, and transport. Qualitatively, geometric modeling offers a basis for molecular visualization, which is crucial for the understanding of molecular structure and interactions. Quantitatively, geometric modeling bridges the gap between molecular information, such as that from X-ray, NMR, and cryo-electron microscopy, and theoretical/mathematical models, such as molecular dynamics, the Poisson-Boltzmann equation, and the Nernst-Planck equation. In this work, we present a family of variational multiscale geometric models for macromolecular systems. Our models are able to combine multiresolution geometric modeling with multiscale electrostatic modeling in a unified variational framework. We discuss a suite of techniques for molecular surface generation, molecular surface meshing, molecular volumetric meshing, and the estimation of Hadwiger's functionals. Emphasis is given to the multiresolution representations of biomolecules and the associated multiscale electrostatic analyses as well as multiresolution curvature characterizations. The resulting fine resolution representations of a biomolecular system enable the detailed analysis of solvent-solute interaction, and ion channel dynamics, whereas our coarse resolution representations highlight the compatibility of protein-ligand bindings and possibility of protein-protein interactions.

Ionic Interactions Are Everywhere

Ionic solutions are dominated by interactions because they must be electrically neutral, but classical theory assumes no interactions. Biological solutions are rather like seawater, concentrated enough so that the diameter of ions also produces important interactions. In my view, the theory of complex fluids is needed to deal with the interacting reality of biological solutions.

Quantum Dynamics in Continuum for Proton Transport I: Basic Formulation

Proton transport is one of the most important and interesting phenomena in living cells. The present work proposes a multiscale/multiphysics model for the understanding of the molecular mechanism of proton transport in transmembrane proteins. We describe proton dynamics quantum mechanically via a density functional approach while implicitly model other solvent ions as a dielectric continuum to reduce the number of degrees of freedom. The densities of all other ions in the solvent are assumed to obey the Boltzmann distribution. The impact of protein molecular structure and its charge polarization on the proton transport is considered explicitly at the atomic level. We formulate a total free energy functional to put proton kinetic and potential energies as well as electrostatic energy of all ions on an equal footing. The variational principle is employed to derive nonlinear governing equations for the proton transport system. Generalized Poisson-Boltzmann equation and Kohn-Sham equation are obtained from the variational framework. Theoretical formulations for the proton density and proton conductance are constructed based on fundamental principles. The molecular surface of the channel protein is utilized to split the discrete protein domain and the continuum solvent domain, and facilitate the multiscale discrete/continuum/quantum descriptions. A number of mathematical algorithms, including the Dirichlet to Neumann mapping, matched interface and boundary method, Gummel iteration, and Krylov space techniques are utilized to implement the proposed model in a computationally efficient manner. The Gramicidin A (GA) channel is used to demonstrate the performance of the proposed proton transport model and validate the efficiency of proposed mathematical algorithms. The electrostatic characteristics of the GA channel is analyzed with a wide range of model parameters. The proton conductances are studied over a number of applied voltages and reference concentrations. A comparison with experimental data verifies the present model predictions and validates the proposed model.

High-order fractional partial differential equation transform for molecular surface construction

Fractional derivative or fractional calculus plays a significant role in theoretical modeling of scientific and engineering problems. However, only relatively low order fractional derivatives are used at present. In general, it is not obvious what role a high fractional derivative can play and how to make use of arbitrarily high-order fractional derivatives. This work introduces arbitrarily high-order fractional partial differential equations (PDEs) to describe fractional hyperdiffusions. The fractional PDEs are constructed via fractional variational principle. A fast fractional Fourier transform (FFFT) is proposed to numerically integrate the high-order fractional PDEs so as to avoid stringent stability constraints in solving high-order evolution PDEs. The proposed high-order fractional PDEs are applied to the surface generation of proteins. We first validate the proposed method with a variety of test examples in two and three-dimensional settings. The impact of high-order fractional derivatives to surface analysis is examined. We also construct fractional PDE transform based on arbitrarily high-order fractional PDEs. We demonstrate that the use of arbitrarily high-order derivatives gives rise to time-frequency localization, the control of the spectral distribution, and the regulation of the spatial resolution in the fractional PDE transform. Consequently, the fractional PDE transform enables the mode decomposition of images, signals, and surfaces. The effect of the propagation time on the quality of resulting molecular surfaces is also studied. Computational efficiency of the present surface generation method is compared with the MSMS approach in Cartesian representation. We further validate the present method by examining some benchmark indicators of macromolecular surfaces, i.e., surface area, surface enclosed volume, surface electrostatic potential and solvation free energy. Extensive numerical experiments and comparison with an established surface model indicate that the proposed high-order fractional PDEs are robust, stable and efficient for biomolecular surface generation.

High-order fractional partial differential equation transform for molecular surface construction

Fractional derivative or fractional calculus plays a significant role in theoretical modeling of scientific and engineering problems. However, only relatively low order fractional derivatives are used at present. In general, it is not obvious what role a high fractional derivative can play and how to make use of arbitrarily high-order fractional derivatives. This work introduces arbitrarily high-order fractional partial differential equations (PDEs) to describe fractional hyperdiffusions. The fractional PDEs are constructed via fractional variational principle. A fast fractional Fourier transform (FFFT) is proposed to numerically integrate the high-order fractional PDEs so as to avoid stringent stability constraints in solving high-order evolution PDEs. The proposed high-order fractional PDEs are applied to the surface generation of proteins. We first validate the proposed method with a variety of test examples in two and three-dimensional settings. The impact of high-order fractional derivatives to surface analysis is examined. We also construct fractional PDE transform based on arbitrarily high-order fractional PDEs. We demonstrate that the use of arbitrarily high-order derivatives gives rise to time-frequency localization, the control of the spectral distribution, and the regulation of the spatial resolution in the fractional PDE transform. Consequently, the fractional PDE transform enables the mode decomposition of images, signals, and surfaces. The effect of the propagation time on the quality of resulting molecular surfaces is also studied. Computational efficiency of the present surface generation method is compared with the MSMS approach in Cartesian representation. We further validate the present method by examining some benchmark indicators of macromolecular surfaces, i.e., surface area, surface enclosed volume, surface electrostatic potential and solvation free energy. Extensive numerical experiments and comparison with an established surface model indicate that the proposed high-order fractional PDEs are robust, stable and efficient for biomolecular surface generation.

Geometric modeling of subcellular structures, organelles, and multiprotein complexes

Recently, the structure, function, stability, and dynamics of subcellular structures, organelles, and multiprotein complexes have emerged as a leading interest in structural biology. Geometric modeling not only provides visualizations of shapes for large biomolecular complexes but also fills the gap between structural information and theoretical modeling, and enables the understanding of function, stability, and dynamics. This paper introduces a suite of computational tools for volumetric data processing, information extraction, surface mesh rendering, geometric measurement, and curvature estimation of biomolecular complexes. Particular emphasis is given to the modeling of cryo-electron microscopy data. Lagrangian-triangle meshes are employed for the surface presentation. On the basis of this representation, algorithms are developed for surface area and surface-enclosed volume calculation, and curvature estimation. Methods for volumetric meshing have also been presented. Because the technological development in computer science and mathematics has led to multiple choices at each stage of the geometric modeling, we discuss the rationales in the design and selection of various algorithms. Analytical models are designed to test the computational accuracy and convergence of proposed algorithms. Finally, we select a set of six cryo-electron microscopy data representing typical subcellular complexes to demonstrate the efficacy of the proposed algorithms in handling biomolecular surfaces and explore their capability of geometric characterization of binding targets. This paper offers a comprehensive protocol for the geometric modeling of subcellular structures, organelles, and multiprotein complexes.