A DNN Biophysics Model with Topological and Electrostatic Features.

A DNN Biophysics Model with Topological and Electrostatic Features - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice J Chem Inf Model . 2026 Mar 23;66(7):3660–3676. doi: 10.1021/acs.jcim.5c03104 Search in PMC Search in PubMed View in NLM Catalog Add to search A DNN Biophysics Model with Topological and Electrostatic Features Elyssa Sliheet Elyssa Sliheet † Department of Mathematics, Southern Methodist University, Dallas, Texas 75275, United States Find articles by Elyssa Sliheet † , Md Abu Talha Md Abu Talha † Department of Mathematics, Southern Methodist University, Dallas, Texas 75275, United States Find articles by Md Abu Talha † , Weihua Geng Weihua Geng † Department of Mathematics, Southern Methodist University, Dallas, Texas 75275, United States Find articles by Weihua Geng †, * Author information Article notes Copyright and License information † Department of Mathematics, Southern Methodist University, Dallas, Texas 75275, United States * Email: [email protected] . Received 2025 Dec 15; Accepted 2026 Mar 3; Revised 2026 Mar 2; Collection date 2026 Apr 13. © 2026 The Authors. Published by American Chemical Society This article is licensed under CC-BY-NC-ND 4.0 PMC Copyright notice PMCID: PMC13080996  PMID: 41870490 Abstract In this project, we present a deep neural network (DNN)-based biophysics model that uses multiscale and uniform topological and electrostatic features to predict protein properties, such as Coulomb energies or solvation energies. The topological features are generated using element-specific persistent homology (ESPH) on a selection of heavy or carbon atoms. The electrostatic features are generated using a novel Cartesian treecode, which adds underlying electrostatic interactions to further improve the model prediction. These features are uniform in number for proteins of varying sizes; therefore, the widely available protein structure databases can be used to train the network. These features are also multiscale, allowing users to balance resolution and computational cost. The optimal model trained on more than 17,000 proteins for predicting Coulomb energy achieves an MSE of approximately 0.024, MAPE of 0.073, and R 2 of 0.976. Meanwhile, the optimal model trained on more than 4000 proteins for predicting solvation energy achieves an MSE of approximately 0.064, MAPE of 0.081, and R 2 of 0.926, showing the efficiency and fidelity of these features in representing the protein structure and force field. The feature generation algorithms also have the potential to serve as general tools for assisting machine learning-based prediction of protein properties and functions. Introduction One of the overarching themes of biology is that structure determines function, which becomes more significant when we explore the biological world at the molecular level. With the advancement and availability of structure determination techniques such as X-ray crystallography, , NMR spectroscopy, cryo-electron microscopy, etc., there are rapidly growing 3D structures of protein, nucleic acids, and complex assemblies, etc. on publicly accessible repository such as Protein Data Bank (PDB), GenBank, etc. These structure determination techniques have strengths and limitations that make them suitable for different types of research studies. For example, X-ray crystallography, due to the requirement of crystallization, is best for static high-resolution structures. NMR is valuable for studying dynamics in solution, particularly for studies of smaller proteins and complexes. Cryo-electron microscopy, which has become more and more popular recently, excels in visualizing large, complex, and flexible biomolecular systems. In addition to experiments, computational approaches are also available for the prediction of protein structures and properties. The main stream of Machine Learning (ML)-based protein prediction focuses on using amino acid sequences to predict protein structure or structure-related properties, resulting in the trending protein language models (PLMs), , which is a transformer-based model inspired by Natural Language Processing (NLP). PLM has a wide range of applications such as protein p K a prediction, protein structure prediction, protein design, protein evolution, etc. Recently, the appearance of AlphaFold developed by Google DeepMind revolutionarily reshaped the structures’ determination techniques. Both PLM and AlphaFold take protein sequences as inputs to predict structures or embeddings that keep important features. However, the development of ML methods in protein property predictions is hindered by the difficulty of representing protein structures and force fields using reliable features. To this end, the structure data, which is largely variant from protein to protein, needs to be converted to a uniform size of features to be loaded into the ML models. Successful cases are the use of graph convolutional network (GCN) for predicting protein functions by features extracted from a protein language model and protein structure, protein–DNA binding specificity prediction using a wide range of features from geometry, , protein–ligand binding affinity prediction using topological features, , molecular structure representation using persistent diagram, etc. However, the important electrostatic interactions are often ignored in the formation of the protein features because of the challenges from the long-range and pairwise nature of the electrostatic interactions. In this paper, we provide a novel approach that utilizes the far-field expansion in the Cartesian treecode algorithm to generate electrostatic features. To this end, the charges at the atomic centers are represented by multipole moments at the cluster centers, which inherit the key idea of the treecode algorithm: replacing pairwise particle–particle interactions with particle–cluster interaction. By combining electrostatic features with the topological features generated using a modification of the persistent homology approach used by Cang and Wei, we represent the protein structures and the associated force field in a uniform and multiscale fashion. With these features, we train a deep neural network (DNN) to predict the protein Coulombic energy and electrostatic solvation energy efficiently and accurately. Here, the Coulomb energy and solvation energy used as labels are generated by pairwise Coulombic interactions and the surrogate model for solving the Poisson–Boltzmann (PB) model using a second-order matched interface and boundary (MIB) method-based PB solver. Apart from the recent research using Physics Informed Neural Network (PINN) to solve the Poisson–Boltzmann model, , the present work is fundamentally different as a data-driven method with uniform and multiscale features. There exist machine learning methods for the prediction of solvation energy under the framework of the PB model. For examples, Huang et al. developed a neural-network formulation for solving the dielectric-boundary PB equation by minimizing a variational electrostatic free-energy functional, enabling direct prediction of PB electrostatic solvation energies while preserving the physical structure of the continuum model; Röcken et al. proposed the ReSolv framework trained using experimental hydration free energies and ab initio data, achieving near experimental accuracy while accelerating solvation energy prediction; Ferraz-Caetano et al. introduced an explainable supervised learning approach based on ensemble regression models that predicts solvation energies from physicochemical descriptors; and Wu and Luo presented SPACIER, a large-scale benchmark dataset of solvation energies grounded computed by using AMBER-PBSA. The major advantage of the proposed method is at the feature generation independent of the PB framework; thus, in addition to this particular applications in predicting Coulombic energy and PB-based electrostatic solvation energy, the feature generation algorithm developed in this work potentially has broader applications to many other ML models, as far as protein structures and force fields are the inputs represented uniformly across proteins of varying sizes. This article is organized as follows. Following this introduction, we present the methods section, which includes both the modeling framework and numerical algorithms, followed by the simulation results. The manuscript ends with a section of the conclusion. Methods The objective of this project is to design uniform and multiscale features that represent protein structure and electrostatics and then to use ML models, such as DNN, to evaluate the effectiveness of these features. The models are trained using labels such as solvation energy computed from surrogate models, e.g., the numerically solved Poisson–Boltzmann model. We hereby introduce the models and algorithms involved in this work. To study protein structure–property relationships, we focus on the atomic description and the electrostatic interactions. The protein properties considered in this paper are Coulomb energy E coul and electrostatic solvation energy E solv . For a protein, the Coulomb energy E coul is the energy required to move atomic charges from infinity to their present positions (i.e., atomic centers specified in the Protein Data Bank (PDB) file). The electrostatic solvation energy E solv is the electrostatic free energy required to solvate the protein in a solvent such as water. Formal quantitative definitions of these two properties are provided through mathematical models. For the protein–solvent interaction, we use the implicit solvent model in which water is treated as a continuum in order to reduce computational cost. From a comparative perspective, we introduce the two most popular implicit solvation models, namely, the Poisson–Boltzmann model and the Generalized Born’s model, which are used to generate labels or features. We then introduce persistent homology and describe how it is used to construct topological features. Finally, we present the Cartesian treecode algorithm, and its modification into a point-multipole representation, which forms the basis for generating electrostatic features. In this project, we adopt a combination of electrostatic and topological features for the following reasons. Electrostatic features provide a uniform and multiscale representation of the atomic charges, accounting for both charge quantities and their locations. Topological features, on the other hand, encode information derived from the atomic centers of selected atoms such as alpha-carbon atoms or heavy atoms. These features reveal the intrinsic topological invariants, which are not directly represented by electrostatic features. Together, these two types of features efficiently represent charge distributions, their locations, and the underlying topological information. The Poisson–Boltzmann Model Figure depicts the popular implicit solvent models. In Figure a, a protein is represented by a collection of N c spherical atoms with centered partial charges. The molecular surface Γ (also known as the solvent excluded surface , ) is defined by the trace of a water molecule represented by a red sphere rolling on contacting with the protein atoms. The Poisson–Boltzmann model is shown in Figure b, where the molecular surface Γ divides the entire computational domain Ω into the protein domain Ω 1 with dielectric constant ε 1 and atomic charges q k located at x k , k = 1: N c , and the solvent domain Ω 2 with dielectric constant ε 2 and dissolved mobile ions. Assuming a Boltzmann distribution for the ion concentration, and considering the case of two ion species with equal and opposite charges (e.g., Na + , Cl – ), in the limit of a weak electrostatic potential one obtains the linearized Poisson–Boltzmann (PB) model, ε 1 ∇ 2 ϕ 1 ( x ) = − 4 π ∑ k = 1 N c q k δ ( x − x k ) , ⁣ x ∈ Ω 1 1a ε 2 ∇ 2 ϕ 2 ( x ) = κ 2 ϕ 2 ( x ) , ⁣ x ∈ Ω 2 1b ϕ 1 ( x ) = ϕ 2 ( x ) , ⁣ ε 1 ∂ ϕ 1 ∂ n ( x ) = ε 2 ∂ ϕ 2 ∂ n ( x ) , ⁣ x ∈ Γ 1c where κ is the inverse Debye length measuring the salt concentration and the potential satisfies a zero far-field boundary condition lim | x | → ∞ ⁡ ϕ ( x ) = 0 . 1. Open in a new tab Implicit solvent models. (a) The molecular surface Γ: the trace of the solvent probe (shown in red) when it is rolled on contacting the spherical atoms of the protein (shown in green with a centered charge). (b) The Poisson–Boltzmann model: two domains Ω 1 (dielectric constant ϵ 1 with partial charges as weighted summation of delta functions) and Ω 2 (dielectric constant ϵ 2 with mobile ions modeled by the Boltzmann distribution) are separated by the molecular surface Γ. (3) The generalized Born model: the protein is represented as a collection of N c spherical atoms centered at x k with charge q k and Born radius R k (only the k th atom is shown). The PB model governs the electrostatic potential ϕ over the entire space. Theoretically after ϕ is obtained, its gradient will produce an electrostatic field while its integral will generate potential energy. However, there are many challenging issues on properly obtaining the field and energy (e.g., definition of field on molecular surface Γ , ). Our attention for this project is on the energy side, as described below. The electrostatic free energy is given as E f r e e = 1 2 ∫ Ω ρ ( x ) ϕ ( x ) d x = 1 2 ∑ k = 1 N c q k ϕ ( x k ) = 1 2 ∑ k = 1 N c q k ( ϕ r e a c ( x k ) + ϕ c o u l ( x k ) ) = E s o l v + E c o u l 2 where ρ ( x ) = ∑ k = 1 N c q k δ ( x − x k ) is the charge density as a sum of partial charges weighted delta functions and the E solv = 1 2 ∑ k = 1 N c q k ϕ r e a c ( x k ) term is the electrostatic solvation energy. The ϕ reac ( x k ) is the reaction potential at x k as the remaining component when Coulomb potential ϕ c o u l ( x k ) = ∑ j = 1 , j ≠ k N c q j ϵ 1 | x k − x j | is taken away from the electrostatic potential ϕ­( x k ). Solving the PB model numerically by grid-based methods is challenging because (1) the protein is represented by singular point charges, (2) the molecular surface is geometrically complex, (3) the dielectric constant is discontinuous across the surface, and (4) the domain is unbounded. To overcome these numerical difficulties, many finite difference interface methods have been developed. − For the purpose of dealing with arbitrarily shaped dielectric interfaces based on a simple Cartesian grid, a matched interface and boundary (MIB) PB solver ,− has been developed through rigorous treatments of geometrical and charge singularities. Alternatively, boundary element methods (BEMs) for the PB model were developed later − with several inherent advantages, (1) only the molecular surface is discretized rather than the entire solute/solvent volume, (2) the atomic charges are treated analytically, (3) the interface conditions are accurately enforced, (4) the far-field boundary condition is imposed analytically, and (5) the matrix-vector product in iterative methods can be accelerated by fast algorithms. ,− In this paper, we computed the solvation energies using the MIBPB solver since accuracy is of critical importance for the generation of labels ( E solv ) that will be used to train the neural network. Generalized Born’s Model If ion effects are not considered, the solvation energy and the reaction potential can also be efficiently computed by using the generalized Born’s (GB) model, an approximation to Poisson’s model. The GB model has significantly reduced computational cost, as opposed to the PB or Poisson models, whose solutions require solving a 3D partial differential equation. Although we did not use the GB model to produce any numerical results for this project, we describe the GB model for its relationship with the PB model and its potential usage in electrostatic interactions. In fact, the GB model can calculate the electrostatic solvation energy which can be used either as the label or as a feature , in our ML framework. It can also produce the reaction potentials at the atomic center, which can be used to generate the second level of electrostatic features, as shown in Figure to be explained within that content. 4. Open in a new tab A 2D illustration for 3D uniform and multiscale electrostatic features for a protein (purple dashed line) with charges (black dots). Charges q and reaction potential ϕ reac are redistributed as point-multipoles (explosion symbols) using Cartesian treecode or FMM at the centers of the cluster at different levels (level 0: black; level 1: red; and level 2: green). Under the GB model, the total electrostatic solvation energy is E s o l v = 1 2 ( 1 ϵ 2 − 1 ϵ 1 ) ( ∑ i = 1 N q i 2 a i + ∑ j ≠ i N q i q j r i j ) ≈ 1 2 ( 1 ϵ 2 − 1 ϵ 1 ) ∑ i = 1 N ∑ j = 1 N q i q j f i j G B 3 where r ij is the distance between the atomic centers of atoms i and j and f ij is the effective Born radii ( i = j ) or effective interaction distance ( i ≠ j ). This PB-GB relationship can be seen in Figure c, in which the electrostatic solvation energy of the sphere equals that of the molecule, with both having a charge at the atomic center. Some details about the GB model can be found in the Supporting Information . Topological Features This section involves many fundamental concepts in Algebraic Topology. We use italic fonts to emphasize important terms and use bold fonts to indicate their sequential definitions, if there are any. The fundamental task of topological data analysis is to extract topological invariants as intrinsic features of the underlying space. For the prediction of protein properties, we expect the topological invariants such as independent components, rings, cavities, etc. carry useful information which cannot be discovered from geometric observation and measurement. This is realized by the aid of homology, a mathematical framework that assigns algebraic objects (groups) to a topological space to measure their shape, especially their holes in various dimensions. Below, we briefly introduce the simplicial homology and the persistent homology and their associated terms and definitions. Following that, we explain how these homologies and their computations with protein structures are used to generate topological features. Simplicial Homology Topological invariants in a discrete dataset (e.g., the collection of atomic positions of a protein) can be studied using simplicial homology which uses a specific rule (e.g., the choice of proximity parameter ϵ and an associated procedure to determine the edges connecting the vertices) to identify simplicial complexes (will be defined later) from simplexes. Here, the simplex represents the simplest possible polytope in any given dimension like a point, line segment, triangle, tetrahedron, etc. Formally, a k -simplex is the convex hull of k + 1 affinely independent points v 0 , v 1 , ..., v k : σ k = { ∑ i = 0 k λ i v i | λ i ≥ 0 , ∑ λ i = 1 } 4 Each subset of m + 1 of these vertices forms an m -face of the simplex, which is itself an m -simplex. The boundary of a k -simplex [ v 0 , ..., v k ] is: ∂ k [ v 0 , v 1 , ... , v k ] = ∑ i = 0 k ( − 1 ) i [ v 0 , ... , v̂ i , ... , v k ] 5 where ∂ k is the boundary operator from k -chains to ( k – 1)-chains and v̂ i means vertex v i is omitted. Here, k -chain is a formal linear combination of k-dimensional simplexes in a simplicial complex. With the concepts of simplex, face, boundary, chain, etc., we define that a simplicial complex K is a collection of simplices such that: 1. Every face of a simplex in K is also in K . 2. The intersection of any two simplices is either empty or a common face. Thus, each simplicial complex corresponds to a topological space by filling in the simplexes and gluing them together along shared faces. Let C k ( K ) be the group of k -chains, and let · · · → ∂ k + 1 C k ( K ) → ∂ k C k − 1 ( K ) → ∂ k − 1 · · · 6 be the sequence of boundary maps ∂ k , the k -homology group of K is defined as: H k ( K ) = ker ( ∂ k ) i m ( ∂ k + 1 ) 7 where ker­(∂ k ) is the group of k -cycles: Z k = { c ∈ C k |∂ k ( c ) = 0} im­(∂ k +1 ) is the group of k -boundaries: B k = {∂ k +1 ( c )| c ∈ C k +1 } Hence, the k -th homology group is the group of cycle modulo boundaries: H k = Z k / B k 8 Taking H 1 as an example, Z 1 = ker­(∂ 1 ) collects 1-chains (formal sum of edges) without boundary points thus general loops, while B 1 = im­(∂ 2 ) collects boundaries of 2-chains (a collection of triangles) that bound a filled-in region thus trivial holes, and H 1 = Z 1 / B 1 is then the collection of actual holes in the space. Persistent Homology However, the simplicial homology of a complex associated with a point set at a particular proximity parameter ϵ is insufficient to characterize the signal or noise. This calls for the ideas of filtration and persistent homology, , which can identify and connect complexes at different levels of complexity and record the appearance and disappearance of the homology group. To this end, we introduce the concept of filtration, which is a sequence of nested simplicial complexes: ϕ = K 0 ⊆ K 1 ⊆ · · · ⊆ K n = K 9 Each K i is a simplicial complex, and the sequence grows over time or on a scale. Instead of computing just the k -homology group H k ( K ), which is an algebraic object capturing the k -dimensional holes in a topological space represented by a simplicial complex K , we compute H k ( K i ) for each i , the index for time or scale parameters. As i or its corresponding scale parameter grows, the duration, as measured by the change of i or its corresponding scale parameter, of the homology class H k ( K i ) between its appearance and disappearance is called persistence. The importance of persistence is that given a parametrized family of spaces, topological features which persist over a significant parameter ranges are considered as signals, while those with short-lived features are treated as noise. In practice, the persistent homology is computed using the pipeline that involves steps of turning data into a sequence of spaces (a filtration), tracking the birth and death of topological features using algebra, implementing the data using linear algebra tools, and visualizing. The pipeline is realized by many established software, e.g., GUDHI, as our choice. We next introduce the barcode, persistent diagram, and persistent image as approaches to visualize the persistent homology. These approaches are different in format and visualization but are equivalent mathematically. Following the definition, we provide the algorithms to extract topological features from barcodes and persistent image. Barcode, Persistence Diagram, and Persistence Image Persistent homology captures the evolution of topological features throughout a filtration process, which could be encoded in the forms of a barcode, persistent diagram, and persistent image. A barcode for dimension k is the multiset { [ b j , d j ) } j , where each interval corresponds to a k -dimensional homology class with index j , which appears at b j and disappears at d j . In the picture of a barcode, each class is visualized as a bar starting at b j and ending at d j . Thus, the barcode is a graphical representation of persistent homology as a collection of horizontal line segments in a plane whose horizontal axis corresponds to the filtration parameter and whose vertical axis represents an (arbitrary) ordering of homology generators such as H 0 , H 1 , H 2 , etc. Different colors can be used to distinguish different dimensions. A persistent diagram represents each homology class as a point ( b j , d j ) in the birth–death plane, where b j and d j denote the filtration parameters at which the class appears and disappears. The persistent diagram thus is a scatter plot of points, where the distance from the diagonal measures their significance. From their definitions, we can tell barcode and persistence diagram carry an equal amount of information. We, therefore, only use a barcode to generate topological features. The persistent image, which is built on top of the persistent diagram to generate a fixed-size 2D image, is readily used as features or inputs for the convolution neural network. Details about persistent image and its formal mathematical definition can be found in the Supporting Information . Feature Generation from Barcode To make the topological features uniform and physically informed, we chose to use element-specific persistent homology (ESPH) to extract topological features at different levels of complexity. ESPH has been successfully applied to the prediction of protein–ligand binding. The core idea is to select atomic elements that form the point cloud by incorporating physical, chemical, and biological information relevant to the target problem. In this project, since solvation energy and Coulomb energy are closely related to the amino acid chains, we use two collections of point clouds. One consists of all carbon atoms, which form the backbone of amino acids. The other consists of all heavy atoms {C, N, O, S}, whose locations and interconnections are also critical to protein structures. Once the collections of point clouds are determined, the corresponding barcodes, which encode the persistent homology of the point clouds, are constructed. We use the GUDHI software package to generate these barcodes. The resulting barcodes are then processed by the algorithm described below to produce vectors of topological features. Figure shows the barcodes of H 0 , H 1 , and H 2 for protein 1a1e, the C-SRC Tyrosine Kinase, and protein 1pb9, the N -methyl- D aspartate receptor, using heavy atoms. The persistent homology analysis conducted in this study characterizes key biophysical phenomena, including intermolecular interactions, solvation effects, and hydrophobicity. Topological features are extracted from H 1 and H 2 barcodes, where the birth, death, and persistence values are computed based on the spatial distribution of carbon atoms and all heavy atoms within the protein structure, resulting in a total of 12 feature channels. 2. Open in a new tab Barcode of H 0 , H 1 , and H 2 generated by the collection of cloud points using heavy atoms in proteins (a) 1a1e and (b) 1pb9. The significance of 1- and 2-dimensional homology groups lies in the fact that they represent 2D rings and 3D voids, respectively, which persist over a wide range of scales and are both topologically and geometrically meaningful. In contrast, the 0-dimensional homology group represents connected components, whose number decreases from the total number of points to 1 and therefore carries a limited persistent signal. Higher-dimensional homology groups, although mathematically defined, are typically short-lived in persistence, behave as combinatorial noise, and lack clear geometric meaning. We use both the collection of heavy atoms and the collection of carbon atoms to capture critical structural and chemical information on proteins. Features derived from all heavy atoms represent geometric interactions between specific types of atoms in the protein, whereas those computed from carbon atoms primarily reflect hydrophobic interactions and indirectly capture solvation effects. The analysis is performed over a distance scale of 0 to 50 Å, which is discretized into uniform bins of 0.25 Å to ensure fine-grained resolution of topological patterns. The topological information encoded in the barcodes is then converted into features using the following approach. We define the collection of barcodes as B ( α , C , D ) with the following specifications. α: atom labels (a specific way to obtain a collection of atoms, i.e., all heavy atoms in protein, heavy atoms close to the ligand, atoms near the mutated residue, etc.) C : type of a simplicial complex (i.e., Ribs or Cech) D : dimension (i.e., H 1 , H 2 , etc.) Using the collection, the structured vectors V b , V d , and V p can be constructed to, respectively, describe the birth, death, and persistent patterns of the barcodes in various spatial dimensions. Practically, the filtration interval [0, L ] is divided into n equal length subintervals, and the patterns are characterized on each subinterval. The description vectors using their i th component for 1 ≤ i < n are defined as: V i b = ∥ { ( b j , d j ) ∈ B ( α , C , D ) | ( i − 1 ) L / n ≤ b j ≤ i L / n } ∥ , 1 ≤ i < n , V i d = ∥ { ( b j , d j ) ∈ B ( α , C , D ) | ( i − 1 ) L / n ≤ d j ≤ i L / n } ∥ , 1 ≤ i < n , V i p = ∥ { ( b j , d j ) ∈ B ( α , C , D ) | ( i − 1 ) L / n ≥ b j , i L / n ≤ d j } ∥ , 1 ≤ i < n , 10 These vectors can be viewed as (1D) images whose pixel value is ∥ ·∥, the cardinality of the set. Each pixel is specialized by i , α, C , D as indices of corresponding sets. In this project, α ∈ {all carbon atoms in the protein, all heavy atoms in the protein}, C = C e c h , D ∈ { H 1 , H 2 } , resulting in 12 vectors with the choice of b , d , p . This barcode-based approach is comparable to the method of persistence imaging, which can be directly integrated into standard machine learning pipelines (e.g., SVMs, neural networks, etc.) at the cost of more preprocessing such as smoothing and discretization as explained previously. Feature Generation from Persistent Image The topological features can also be constructed using a persistence image alternatively. To this end, we first extract persistence diagrams from protein structures and convert them to stable, fixed-size image representations. The persistence diagrams are obtained separately from two types of point clouds: the collection of all carbon atoms and the collection of all heavy atoms. These two approaches are also used in the barcode-based feature construction, as detailed previously. For each approach, we consider topological features in homology groups H 1 and H 2 . Each persistence diagram is filtered by using physically meaningful thresholds such that birth and persistence values are restricted to maximum bounds of 12.5 Å. The resulting persistence domain is discretized using a resolution of 0.25 Å, which produces a uniform image grid D = [0, b thr ] × [0, p thr ] with b thr = p thr = 50 for all proteins. To emphasize the contribution of long-lived topological features, each point ( b i , p i ) ∈ D is assigned a smooth exponential weight ω ( b i , p i ) = e p i − 1 11 This exponential weighting enhances features with larger persistence while assigning near-zero weight to short-lived topological noise, thereby improving the contrast between significant and insignificant structures in the resulting persistence images. In this project, persistence images are generated using the Gaussian kernel with a bandwidth of 0.4 Å, which is implemented by using the GUDHI library. For each protein, four persistence images are produced: carbon atoms only for H 1 and H 2 and all heavy atoms for H 1 and H 2 . These four images are stacked to form a four-channel topological feature tensor, X ∈ R b t h r × p t h r × 4 . Finally, the resulting persistence image tensor is stored in the NumPy format and used as the input topological representation for downstream deep learning models. This approach provides a stable, compact, and fully vectorized encoding of protein topology and thus is well suited for ML-based prediction. In Figure , the barcodes of H 1 and H 2 using carbon atoms and using all heavy atoms for proteins 1a0t and 1a4q are provided for comparison purposes with the persistent diagram and persistent image. The persistence diagrams of the two proteins 1a0t and 1a4q show that most of the topological features are short-lived and concentrated near the diagonal, while only a small number of them persist. The corresponding persistence images highlight these persistent structures as bright and localized regions, with background noise effectively suppressed. The observed differences between persistence images from the cloud of points of carbon atoms and those from all heavy atoms show that topological invariants are influenced by the underlying chemical composition. 3. Open in a new tab Barcode (top), persistence diagrams (middle), and persistence image (bottom) for proteins 1a0t (a) and 1a4q (b). For each protein, from left to right are H 1 and H 2 for carbon atoms and H 1 and H 2 for all heavy atoms. Electrostatic Features Introduction In most molecular simulations, including Monte Carlo simulation, Brownian dynamics, molecular dynamics, etc., the electrostatic interactions are characterized by the partial charges assigned at the atomic centers. These partial charges are assigned according to the force fields, which are a set of mathematical equations and parameters used to calculate the potential energy of a molecule based on its atomic structure and interactions. The force fields are determined by experiment or quantum chemistry. The algorithm for obtaining the multiscale, physics-informed, uniform electrostatic features can be explained using Figure . The electrostatic profile of the protein is determined by the partial charges q n ( r n ) or reaction potential ϕ reac ( r n ) at the atomic centers r n ’s as illustrated by the black dots for n = 1, 2, ..., N c with N c as the number of atoms of the protein. Partial charges are the source of electrostatic interactions, while reaction potentials are the outcomes of electrostatic interactions for the solvated protein, which includes the important solute–solvent interactions. Considering the variation of N c , the number of atoms, from protein to protein, it is difficult to directly use q n ( r n ) and/or ϕ reac ( r n ) as features for machine learning. Our strategy to alleviate this pain is to use point-multipoles, as illustrated by the explosion symbols at the cluster centers, to carry redistributed charges or reaction potentials in terms of moments of the multipole expansion. The number of cluster centers is uniform for proteins with various numbers of atoms. It is also multiscale as being determined by the combination of levels of the tree L and number of terms of the multipole expansion p . The user specifies L and p as needed. In Figure , the N d point multipoles (explosion symbols) are organized in a hierarchical way such that different levels of details are given by different numbers of symbols at different levels. For the m th multipole, m = 1, ..., N d , using a 3D index k = ( k 1 , k 2 , k 3 ), | k | = k 1 + k 2 + k 3 , in a cluster c with cluster center x c , the moments are defined as M c k = ∑ x j ∈ c q j ( x j − x c ) k 12 Note when | k | = 0, k = (0, 0, 0); when | k | = 1, k ∈ {(1, 0, 0), (0, 1, 0), (1, 0, 0)}; when | k | = 2, k ∈ {(1, 1, 0), (1, 0, 1), (0, 1, 1), (2, 0, 0), (0, 2, 0), (0, 0, 2)}, etc. If the Cartesian treecode algorithm is used, the computational cost is O ( N c log N c ) for obtaining all the moments as detailed in where N c is number of the particles/charges and log N c is the level of trees. In fact, the computational cost can be as low as O ( N c ) using the strategies in such that the moments at the finest cluster are computed first and a M2M (moments to moments) transformation can be used to efficiently compute moments at any desired level from the bottom to the top. Our simulations so far showed that the O ( N c log N c ) cost is sufficiently fast for moments computation since N c is within the scale of tens of thousands even for large proteins. These moments are intrinsic properties of the cluster and thus can serve as features for the protein, which carries simplified but important structural and charge information. The number of features up to level L cluster (the largest box containing all particles is named the level 0 cluster and each subdivision increases the level by 1) is N f ( p , L ) = N p ( 1 + 8 + 8 2 + · · · + 8 L ) = N p 8 L + 1 − 1 7 13 where N p is the number of terms in multipole expansion with N p = 1, 4, 10, 20, 35, 56··· as the sequence of tetrahedral numbers when p th order multipole is used. Implementation Details We provide a step-by-step procedure by which the electrostatic features are generated. Step 1. The protein structure (e.g., the PDB file) is obtained from publicly available repo, e.g., the protein data bank. The protein is represented by its atomic locations and partial charges: q i (r i ) for i = 1, ..., N c , e.g., the PQR file using the PDB 2PQR software. Step 2. Two parameters p and L are given, where p is the order in the Cartesian Taylor expansion and L is the number of levels. With given p , we have the number of terms given as N p = ( p + 1)­( p + 2)­( p + 3)/6. The number of features up to level L cluster is N f ( p , L ) = p ( p + 1 ) ( p + 2 ) ( 8 L + 1 − 1 ) 42 14 These N f numbers are ordered by level from 0 to L and the coordinates of cluster centers in each level. This determines the dimension of our final feature vector F ∈ R N f . The following table gives how N f varies with the change of p and L : p\L 0 1 2 3 4 0 1 9 73 585 4681 1 4 36 292 2340 18,724 2 10 90 730 5850 46,810 3 20 180 1460 11,700 93,620 4 35 315 2555 20,475 163,835 Open in a new tab Step 3. The Cartesian treecode will be called and the moments for the multipoles will be output as the electrostatic features. These features can be used as part of the inputs for ML algorithms, which will be shown next. The Machine Learning Model The main goal of this project is to discover the hidden and useful information embedded in the protein structural data from the protein database in a simple and abstract way. The protein structure data and force field include bonded and nonbonded interactions with which the molecular simulation can be performed. We categorize these information into electrostatic, topological, algebraic, and geometric features. These four aspects are all important to the machine learning model. In this project, we put our focus on electrostatic and topological features. The DNN based ML model is shown in Figure , which uses the available protein structural data, force field such as AMBER, CHARMM, AMOEBA, etc. and known protein properties’ repositories such as PDBbind Database, Protein p K a Database, etc. to mathematically generate features to train a learned model and then use the learned model to predict unknown properties for protein with available structures. In this project, we focus on the effectiveness of topological and electrostatic features using computed Coulomb energies and solvation energies as labels. 5. Open in a new tab The DNN-based machine learning model using protein structures, force field, and known protein properties to mathematically generate electrostatic, topological, algebraic, and geometric features to train a learned model and then use the learned model to predict unknown properties for protein with available structures. This project focuses on electrostatic and topological features. Mathematically, the surrogate models (Coulomb energies computed by treecode or solvation energies from the numerically solved PB model) E = f ( X , Q ) 15 are approximated by the machine learning models Ê = f̂ ( F e , F t , Θ ) 16 where E and Ê are the energies, X and Q are the positions and charges of protein atoms, Θ is the parameters in the ML model to be trained by data, and F e = F e ( X , Q ) and F t = F t ( X ) are electrostatic features and topological features and their generators, which are nontrivial mappings from protein structure and electrostatics to features. Specification of the ML Model The model developed in this study is to predict Coulomb and solvation energies using combined topological and electrostatic features. All target labels were standardized to zero mean and unit variance, and predictions were inversely transformed to recover physical units (kcal/mol) for evaluation. As illustrated in Figure , the proposed architecture consists of two branches. The first branch, shown in blue, is a 1D convolutional neural network (CNN) designed to process the topological features. These features are passed through three consecutive convolution–pooling blocks (with 128, 64, and 64 filters; kernel size = 3; and activation functions tanh, ReLU, and tanh, respectively), followed by a fully connected layer with 64 units. The electrostatic features are processed through two fully connected hidden layers with 128 and 64 units, respectively (shown in orange). After independent processing in these two branches, their outputs are concatenated and passed through three additional fully connected layers with 32, 16, and 8 units (shown in purple). The final dense layer then produces the regression output for the predicted energy (shown in red). To improve generalization, dropout, l 2 regularization, and batch normalization were applied. The model was trained using the Adam optimizer (learning rate of 1.0 × 10 –4 ) with mean squared error as the loss function. The dataset was randomly divided into 80% training and 20% testing subsets. Model selection was performed using 5-fold cross-validation on the training set, and the best-validated model was subsequently evaluated on the held-out test set. The model architecture and hyperparameters were selected based on empirical considerations and adopted from previously published studies. The complete hyperparameter configuration is provided in Table S1 of the Supporting Information 6. Open in a new tab The architecture of the deep neural network with two branches taking in topological features and electrostatic features. The outputs of the two branches are concatenated into one deep neural network. Data The datasets used in this study are sourced from the PDBbind database, a widely used benchmark for protein–ligand complexes. We categorize the data into two datasets as described below. For all structures, atomic van der Waals radii and partial charges are assigned using the AMBER ff14SB force field through PDB 2PQR software. Dataset 1. This dataset is derived from the refined set of the PDBbind v2018 release, which contains 4,463 protein structures. Outliers in the Coulombic energy E coul and solvation free energy E solv are removed using the interquartile range (IQR) method, in which samples falling outside the interval [ Q 1 – 1.5 × IQR, Q 3 + 1.5 × IQR] are excluded. For E coul , the filtering thresholds are [−817.87, 170.83], where Q 1 = – 447.11 kcal/mol, Q 3 = – 199.94 kcal/mol, and IQR = 247.17 kcal/mol. For E solv , the thresholds are [−10, 691.76, 2, 210.83], where Q 1 = −5, 853.29 kcal/mol, Q 3 = −2, 627.64 kcal/mol, and IQR = 3, 225.65 kcal/mol. After this filtering step, 4,043 and 4,088 protein structures are retained for E coul and E solv , respectively, and used for model training and evaluation. All E solv values are computed using the second-order accurate Matched Interface and Boundary Poisson–Boltzmann (MIBPB) solver with a mesh spacing of h = 0.2 Å. Dataset 2. This dataset is constructed from the PDBbind v2020 release, comprising a general set of 14,127 complexes and a refined set of 5,316 complexes. Coulombic energies E coul are computed for all proteins using a Cartesian treecode algorithm. To ensure data quality, outliers are identified and removed using the IQR method, applying thresholds of [−1, 164.23, 338.04] kcal/mol, where Q 1 = −600.88 kcal/mol, Q 3 = −225.31 kcal/mol, and IQR = 375.57 kcal/mol. Following this filtering step, a total of 17,007 protein structures are retained for subsequent model training and evaluation. Figure illustrates the distributions of Coulombic and solvation energies across the two datasets. The Coulombic energies in both datasets are concentrated within a relatively narrow negative range, whereas the solvation energies exhibit a substantially broader distribution with greater dispersion. This contrast reflects the fundamentally different physical scales and variability associated with Coulombic interactions and solvation effects. 7. Open in a new tab Box-and-whisker plots showing the distributions of Coulombic energy E coul for Dataset 1 in (a) and Dataset 2 in (b), and solvation energy E solv for Dataset 1 in (c). Jittered points represent individual protein samples. Research Purposes and Metrics We provide simulation results to evaluate the performance of the proposed model. These results are designed to address the following three research purposes: Purpose 1: Training on a larger dataset leads to improved predictive performance. Purpose 2: When electrostatic features are incorporated, model performance improves as the treecode level L , Taylor expansion order p , and corresponding number of features N f increase. Purpose 3: Combining electrostatic and topological features outperforms models that use only one type of feature. Metrics: We use three common metrics to measure the performance of the model, as defined below. In all these metrics, we use y i to denote actual value, ŷ i for predicted value, y̅ for mean of actual values, and y ^ ® for mean of predicted values. MAPE: The mean average percentage errors, which measures prediction accuracy as a percentage. This metric is not suitable if the actual values contain zero or are close to zero. M A P E = 1 n ∑ i = 1 n | y i − ŷ i y i | × 100 % 17 MSE: The mean-squared error measures the average squared difference between predicted and actual values. M S E = 1 n ∑ i = 1 n ( y i − ŷ i ) 2 18 R 2 , the quantity indicating the proportion of variance in the dependent variable that is predictable from the independent variables. R 2 = 1 − S S r e s S S t o t 19 where S S t o t = ∑ i = 1 n ( y i − y̅ ) 2 measures the total variance in the observed data and S S r e s = ∑ i = 1 n ( y i − ŷ i ) 2 measures the variance not explained by the model. Results We next provide simulation results using different datasets and models. Following that, additional comparisons with baseline models and different grouping strategies are provided to justify the learned model. Model Performance In reporting and explaining the performance of each combination of model and dataset, we focus on achieving the three aforementioned purposes. In describing the PB and GB models, we use the popular e c /Å units for potential and e c /Å units for energy. By multiplying the conversion factor 332.0716, the units of the potential and energy ( E coul and E solv ) become kcal/mol/ e c and kcal/mol, respectively, which are the units reported for all results. E coul Predicted Using Electrostatic Features Only on Dataset 1 and Dataset 2 We first investigated the effectiveness of electrostatic features in predicting the Coulombic energy E coul . Note that E coul is determined by the locations and charges of the particles, which also determine the moments of the multipoles as electrostatic features. Table reports the E coul predictions obtained using electrostatic features only on Dataset 1 and Dataset 2. The first important pattern observed from the table is that the MSE and R 2 obtained from Dataset 2 are generally better than those from Dataset 1, as the differences in MSE (D2–D1) are almost all negative, while the difference in R 2 (D2–D1) is almost all positive. This pattern validates Purpose 2, indicating the advantage of using a larger dataset, particularly in capturing more overall variance in the target. We also noticed that the differences in MAPE­(D2–D1) are almost all positive, apparently contradicting the trends observed for MSE and R 2 . In fact, this behavior provides additional insight into the data distribution: Dataset 2, but not Dataset 1, contains more samples with smaller target values, as shown in Figure , leading to larger percentage errors. 1. E coul (kcal/mol) Prediction Results Using Electrostatic Features Only on Dataset 1 and Dataset 2 . p L N f ( p , L ) dataset 1 dataset 2 D2–D1 MSE MAPE R 2 MSE MAPE R 2 MSE MAPE R 2 1 0 4 0.584 0.337 0.385 0.682 0.531 0.323 0.098 0.194 –0.062 2 0 10 0.353 0.246 0.628 0.340 0.331 0.663 –0.013 0.085 0.035 3 0 20 0.254 0.207 0.733 0.236 0.284 0.767 –0.017 0.077 0.035 4 0 35 0.325 0.186 0.668 0.249 0.321 0.753 –0.076 0.135 0.085 0 1 9 0.708 0.344 0.254 0.633 0.498 0.372 –0.075 0.154 0.118 1 1 36 0.256 0.205 0.730 0.241 0.273 0.761 –0.015 0.068 0.030 2 1 90 0.219 0.177 0.769 0.156 0.230 0.845 –0.063 0.053 0.076 3 1 180 0.169 0.155 0.822 0.163 0.237 0.838 –0.006 0.083 0.016 4 1 315 0.185 0.173 0.805 0.165 0.234 0.837 –0.020 0.061 0.032 0 2 73 0.403 0.278 0.575 0.322 0.324 0.681 –0.081 0.047 0.105 1 2 292 0.234 0.223 0.754 0.127 0.206 0.874 –0.106 –0.017 0.120 2 2 730 0.204 0.187 0.786 0.133 0.208 0.868 –0.070 0.021 0.082 3 2 1460 0.182 0.163 0.808 0.123 0.195 0.880 –0.059 0.032 0.071 4 2 2555 0.247 0.225 0.740 0.145 0.229 0.857 –0.102 0.004 0.118 0 3 585 0.393 0.239 0.586 0.265 0.279 0.737 –0.128 0.041 0.151 1 3 2340 0.443 0.211 0.533 0.253 0.236 0.752 –0.190 0.024 0.219 2 3 5850 0.411 0.206 0.567 0.217 0.252 0.784 –0.194 0.046 0.218 3 3 11,700 0.472 0.249 0.503 0.209 0.261 0.792 –0.262 0.012 0.290 4 3 20,475 0.381 0.240 0.598 0.223 0.267 0.779 –0.158 0.027 0.180 Open in a new tab a The D2–D1 columns report the differences (Dataset 2–Dataset 1) for the corresponding metrics. From Table , we also observe that in Dataset 1, the MSE decreases and R 2 increases as p and L increase, reaching an optimal value at p = 3 and L = 1 with 180 features. Similarly, in Dataset 2, the MSE decreases and R 2 increases with increasing p and L , reaching an optimal value at p = 3 and L = 2 with 1460 features. This pattern validates Purpose 1, demonstrating the advantage of using more features. However, we should be aware that the optimal choices of p and L are influenced by the size of the training dataset. When the number of features exceeds the number of samples, increasing the number of features may have a negative impact, such as overfitting. We plot the R 2 values in Figure a, where the horizontal axis represents the number of features N f given in the third column of Table . The validation of Purpose 1 is evident, as the blue squares are consistently lower than the red ones. The validation of Purpose 2 is also observed, as markers of the same color exhibit an increasing trend with the number of features until the feature set becomes excessively large. We also plot the best-performing cases for both datasets in Figure b,c, where the reference-predicted values are distributed closely around the ŷ = y line. 8. Open in a new tab (a) R 2 of model performance for models trained on the two datasets against number of features; (b) scatter plot for Dataset 1 ( p = 3, L = 1 with total 180 features); and (c) scatter plot for Dataset 2 ( p = 3, L = 2 with total 1460 features). E coul Predicted Using e-Features and Both Features on Dataset 2 We next investigated the effectiveness of the combined feature representation. This analysis consists of two parts. In the first part, E coul is used as the label on Dataset 2, while in the second part, E solv is used as the label on Dataset 1, taking into account the size of the datasets and the availability of labels. The MSE, MAPE, and R 2 values for the different models are reported in Table . Using topological features alone yields MSE, MAPE, and R 2 values of 0.030, 0.071, and 0.970, respectively, indicating a strong predictive accuracy. Moreover, the difference column clearly shows that incorporating topological features into the electrostatics-only model consistently reduces MSE and MAPE while increasing R 2 . These results partially support Purpose 2. 2. E coul (kcal/mol) Prediction Results Using Electrostatic Features and Combined Features on Dataset 2 . p L N f ( p , L ) e-features only both features Difference MSE MAPE R 2 MSE MAPE R 2 MSE MAPE R 2 1 0 4 0.682 0.531 0.323 0.028 0.076 0.972 –0.654 –0.454 0.649 2 0 10 0.340 0.331 0.663 0.035 0.072 0.966 –0.305 –0.258 0.303 3 0 20 0.236 0.284 0.767 0.031 0.064 0.969 –0.205 –0.220 0.202 4 0 35 0.249 0.321 0.753 0.027 0.082 0.973 –0.222 –0.239 0.221 0 1 9 0.633 0.498 0.372 0.029 0.067 0.971 –0.604 –0.431 0.599 1 1 36 0.241 0.273 0.761 0.030 0.068 0.971 –0.211 –0.205 0.210 2 1 90 0.156 0.230 0.845 0.029 0.063 0.971 –0.127 –0.167 0.126 3 1 180 0.163 0.237 0.838 0.029 0.067 0.971 –0.134 –0.170 0.133 4 1 315 0.165 0.234 0.837 0.024 0.073 0.976 –0.141 –0.162 0.140 0 2 73 0.322 0.324 0.681 0.031 0.069 0.969 –0.291 –0.255 0.289 1 2 292 0.127 0.206 0.874 0.035 0.131 0.965 –0.092 –0.075 0.092 2 2 730 0.133 0.208 0.868 0.033 0.109 0.968 –0.101 –0.099 0.100 3 2 1460 0.123 0.195 0.880 0.033 0.092 0.967 –0.089 –0.103 0.087 4 2 2555 0.145 0.229 0.857 0.026 0.067 0.974 –0.119 –0.162 0.117 0 3 585 0.265 0.279 0.737 0.031 0.071 0.969 –0.234 –0.209 0.232 1 3 2340 0.253 0.236 0.752 0.034 0.077 0.967 –0.219 –0.159 0.215 2 3 5850 0.217 0.252 0.784 0.031 0.075 0.969 –0.186 –0.177 0.185 3 3 11,700 0.209 0.261 0.792 0.029 0.083 0.971 –0.181 –0.178 0.179 4 3 20,475 0.223 0.267 0.779 0.023 0.074 0.978 –0.200 –0.193 0.199 Open in a new tab a The corresponding MSE, MAPE, and R 2 values obtained using Topological Features only are 0.030, 0.071, and 0.970, respectively. If we consider topological features first, it appears that adding electrostatic features does not significantly improve the results. This is because the performance obtained by using topological features is already very strong, leaving limited room for the effectiveness of electrostatic features to manifest. Our next test case shows that when the topological features are relatively less effective, adding electrostatic features can significantly improve performance. Before moving to that case, we briefly examine Figure , which shows scatter plots of predicted values versus true values using topological features, electrostatic features, and combined features. The improvement from (b) to (c) is very evident, while no obvious improvement is observed from (a) to (c). 9. Open in a new tab Scatter plots for the prediction of E coul on Dataset 2: (a) topological features only, (b) e-features only ( p = 3, L = 2 with total 1460 features), and (c) both features ( p = 4, L = 1 with total 315 features). E solv Predicted Using e-Features Only and Both Features on Dataset 1 As mentioned earlier, computing electrostatic solvation energies accurately is time-consuming. Here, we use the MIBPB solver to compute the solvation energies of proteins in Dataset 1, which are then used as labels for model training and testing. A similar simulation using Dataset 1 is conducted, and the results are reported in Table . From the table, we can see that the improvement obtained by using both features, compared with using only topological features or only electrostatic features, is significant. If we focus on the difference column, we see that adding topological features clearly improves the electrostatics-only model by reducing the MSE and MAPE while increasing the R 2 score. In addition, the results obtained using both features are noticeably better than those obtained using topological features alone, which yield MSE = 0.150, MAPE = 0.153, and R 2 = 0.826. These results validate Purpose 3 and demonstrate the advantage of combining topological and electrostatic features. 3. E solv (kcal/mol) Computed Using e-Features and Both Features on Dataset 1; the MSE, MAPE and R 2 for Using Topological Features Only is 0.150, 0.153 and 0.826 Respectively. p L N f ( p , L ) e-features only both features Difference MSE MAPE R 2 MSE MAPE R 2 MSE MAPE R 2 1 0 4 0.345 0.243 0.604 0.063 0.108 0.926 –0.282 –0.136 0.323 2 0 10 0.271 0.256 0.690 0.090 0.094 0.895 –0.180 –0.162 0.206 3 0 20 0.135 0.142 0.845 0.064 0.081 0.926 –0.071 –0.061 0.081 4 0 35 0.121 0.138 0.861 0.065 0.083 0.925 –0.056 –0.055 0.064 0 1 9 0.447 0.318 0.487 0.114 0.169 0.867 –0.333 –0.149 0.380 1 1 36 0.148 0.161 0.830 0.077 0.104 0.911 –0.072 –0.057 0.081 2 1 90 0.095 0.130 0.891 0.084 0.098 0.903 –0.012 –0.032 0.012 3 1 180 0.111 0.138 0.873 0.105 0.103 0.878 –0.006 –0.036 0.005 4 1 315 0.128 0.144 0.853 0.097 0.105 0.887 –0.031 –0.039 0.034 0 2 73 0.214 0.205 0.754 0.092 0.116 0.892 –0.122 –0.089 0.138 1 2 292 0.169 0.156 0.806 0.127 0.123 0.852 –0.041 –0.033 0.045 2 2 730 0.164 0.156 0.812 0.128 0.116 0.851 –0.036 –0.039 0.039 3 2 1460 0.212 0.168 0.757 0.138 0.124 0.839 –0.074 –0.045 0.082 4 2 2555 0.218 0.168 0.750 0.165 0.128 0.808 –0.053 –0.040 0.058 0 3 585 0.336 0.211 0.615 0.148 0.135 0.827 –0.188 –0.076 0.213 1 3 2340 0.399 0.216 0.542 0.218 0.143 0.747 –0.181 –0.073 0.205 2 3 5850 0.390 0.202 0.553 0.171 0.139 0.801 –0.219 –0.063 0.249 3 3 11,700 0.469 0.242 0.461 0.195 0.145 0.773 –0.275 –0.097 0.312 4 3 20,475 0.365 0.210 0.581 0.242 0.156 0.718 –0.123 –0.053 0.137 Open in a new tab The optimized case for the electrostatic-only model occurs at p = 2 and L = 1 with 90 features. Adding topological features further improves the performance, achieving an MSE of 0.064, an MAPE of 0.081, and an R 2 value of 0.926 using only 20 features at p = 3 and L = 0. Figure presents scatter plots for the optimized cases, comparing reference and predicted values for models using topological features only, electrostatic features only, and the combined feature set. The benefits of combining both types of features, as illustrated in (c), can be clearly observed when compared with using only topological features in (a) and only electrostatic features in (b). 10. Open in a new tab Scatter plots for the prediction of E solv on Dataset 1: (a) topological features only, (b) e-features only ( p = 2, L = 1 with total 90 features), and (c) both features ( p = 3, L = 0 with total 20 e-features). Additional Comparison and Justification We here provide some comparison studies between the proposed models with baseline models, surrogate models, and different grouping and persistent homologies’ strategies. In addition, some details about the cross-validation in the training process from a selected case are provided. The limitations of the proposed model and the feature generation algorithms are also discussed. E coul Predicted from Persistence Image on Dataset 2 As explained previously, the topological features can be extracted using either barcodes or persistent images. When only the persistence image-based topological features are used, the model achieves an MSE of 0.037, MAPE of 0.082, and an R 2 value of 0.964, which is comparable to the performance obtained using barcode-based topological features (MSE = 0.030, MAPE = 0.071, and R 2 = 0.970). Both approaches demonstrate excellent predictive accuracy, even in the absence of electrostatic features. When these topological features are combined with electrostatic features, the barcode-based representation yields slightly better performance than the persistent image-based representation, as shown in Figure , in terms of R 2 , MSE, and MAPE. Considering its relatively better performance, lower memory usage, and more efficient computation, we adopt the barcode-based approach throughout the paper. It is worth noting, however, that the persistent image approach may offer advantages in terms of stability and vectorized structure, which will be investigated further in future work. 11. Open in a new tab Prediction performance of E coul on Dataset 2 using combined features: (a) R 2 score, (b) scaled MSE, and (c) scaled MAPE as functions of the number of features, comparing Barcode and PImage representations. Robustness of E coul Prediction to Homology-Aware Data Splitting (Dataset 2) In our machine learning model, we first randomly divide the dataset into training sets and testing sets with proportions of 80% and 20%, respectively. The training set is further divided into five equal subsets for 5-fold cross-validation, from which the model with the best validation performance is selected and subsequently evaluated on the held-out test set. A potential concern is that the PDBbind dataset contains multiple structures corresponding to the same protein or closely related homologues (e.g., alternative conformations, mutants, or orthologous proteins). If such similar proteins are distributed in both training and testing sets, data leakage may artificially inflate the performance. To mitigate this risk, we adopted two complementary strategies. First, we grouped PDB entries according to their UniProt (Universal Protein Resource) accession, ensuring that all structures of the same protein were assigned entirely to either the training or testing set. However, since proteins with different UniProt IDs may still share high sequence similarity, this alone does not eliminate homologous overlap. We therefore applied MMseqs (Many-against-Many sequence searching) clustering at 30% sequence identity and 80% coverage. Chain sequences were clustered, and PDB entries were grouped by connected components such that any entries sharing sequence similarity, including transitive relationships, were placed in the same group. Dataset splitting was then performed at the group level to prevent sequence-similarity leakage. As shown in Figure , similar trends in R 2 , scaled MSE, and scaled MAPE are observed across random, UniProt-based, and MMseqs-based splits. The slightly better performance under random splitting is expected due to less stringent generalization. Comparable results under these methods indicate that the model’s performance is not significantly driven by the memorization of homologous sequences. 12. Open in a new tab Performance comparison for E coul prediction on Dataset 2 under random splitting, UniProt-based grouping, and MMseqs sequence-identity clustering (30% identity, 80% coverage), using combined features. Cross-Validation Performance Now, we investigate the training process that uses 5-fold cross-validation for overfitting checks. In Figure , we report the cross-validation training histories or learning curves for E coul prediction ( p = 4 and L = 1) on Dataset 2. From the mean training and validation loss curves and the corresponding one-standard deviation bands, we can see that the selected model is stable and free of overfitting. Meanwhile, the performance of using both types of features is significantly better than that of using electrostatic features only, as indicated by faster convergence and narrower bands. Similar trends are observed for models predicting other physical quantities across different datasets. 13. Open in a new tab Cross-validation training histories for E coul prediction ( p = 4 and L = 1) on Dataset 2. (a) Electrostatic features only and (b) Combined features. The mean training and validation losses are shown as solid lines, with shaded bands representing standard deviation across folds. We further show in Figure the mean, 95% confidence interval (assuming t -distribution), and one standard deviation band of R 2 in (a) and MSE in (b) from the validation results of the 5-fold cross-validation training procedure across different numbers of electrostatic features when combined features are used. The narrow band and interval widths demonstrate the stability of the model. Note that for E coul , the advantages of using more electrostatic features are not obvious since the contribution of topological features dominates. 14. Open in a new tab Validation performance during training across all cross-validation folds using combined features. Solid lines show mean performance across folds, while shaded regions and error bars indicate 95% confidence intervals and ± one standard deviation, respectively. Comparison with Baseline Models We next present the comparison results between the deep neural network model and other baseline models, including ridge regression, random forest, and gradient boosting. All methods use the same number of electrostatic and topological features, with the number of electrostatic features increasing as shown in Figure . From the comparison, we observe that DNN, random forest, and gradient boosting, as more flexible methods, perform significantly better than ridge regression. The DNN shows slight improvements over random forest and gradient boosting in terms of R 2 and MSE. 15. Open in a new tab Prediction performance of the proposed DNN model compared with baseline methodsridge regression, random forest, and gradient boostingacross varying numbers of input features. CPU Time Comparison and Size-Dependent Performance In order to show the efficiency of the learned model for predicting E solv compared with the MIBPB solver, we report the CPU time and solvation energies E solv obtained from both approaches as shown in Figure . As shown in (a), the predicted E solv values are closely aligned; whereas in (b), the time required by the learned model in red is significantly shorter than the runtime of the MIBPB solver in black, which increases rapidly with molecular size. Also as shown in Figure b, the total machine learning pipeline time in red is further decomposed into time to compute electrostatic features in blue, time to compute topological features in orange, and the time to predict solvation energies using the learned model in green. 16. Open in a new tab Performance comparison between the proposed DNN model and the MIBPB solver. (a) E solv predicted by the learned model versus that is computed by MIBPB solver selected protein structures and (b) computational runtime with respect to protein size (number of atoms), demonstrating the significant efficiency of the DNN model. The MIBPB solver is evaluated with grid size h = 0.5, dielectric constants ϵ 1 = 1, ϵ 2 = 80, and inverse Debye length κ = 0. To further evaluate the robustness of the model with respect to protein size, Figure presents a size-dependent performance analysis. The scatter plot (predicted value vs reference value) in (a) and Box-and-whisker plot (absolute percentage error) in (b) both indicate the model’s uniform performance across protein sizes. 17. Open in a new tab Protein size-dependent performance of the DNN model. (a) Predicted versus reference E solv colored by protein size. (b) Distribution of absolute percentage error across protein size ranges. Model Comparison with Testing Hypothesis Theoretically, the three research purposes as mentioned earlier can be treated as statistical hypotheses and evaluated through testing hypotheses for model comparison. However, this would require pairwise comparisons between models with numerous simulations to generate statistically significant metrics, which is a costly procedure. Here, we present one such case to illustrate the approach. We select a particular set of electrostatic features ( p = 2, L = 1) to compare three models M e , M t , and M e,c , which use electrostatic features only, topological features only, and combined features, respectively. These three models are used to predict E solv on Dataset 2 using 10 independent random 80–20 train–test splits, and the corresponding R 2 and MSE values are recorded. The resulting R 2 and MSE metrics, summarized in Table , are then used to test the following hypotheses: H 0 : M e,t performs as well as M e or M t H a : M e,t performs better than M e or M t 4. Metrics Obtained by Running Models M e , M t , and M e , t for the Prediction of E solv Using 10 Independent Random 80–20 Train–Test Splits on Dataset 2 . metrics S 1 S 2 S 3 S 4 S 5 S 6 S 7 S 8 S 9 S 10 S i ® R 2 ( M e ) 0.876 0.885 0.883 0.853 0.876 0.859 0.879 0.862 0.873 0.885 0.873 R 2 ( M t ) 0.838 0.802 0.816 0.847 0.857 0.803 0.844 0.812 0.852 0.810 0.828 R 2 ( M e , t ) 0.935 0.936 0.907 0.937 0.926 0.887 0.931 0.914 0.921 0.899 0.919 MSE( M e ) 0.120 0.121 0.123 0.147 0.116 0.148 0.125 0.142 0.126 0.111 0.128 MSE( M t ) 0.157 0.214 0.199 0.155 0.133 0.207 0.161 0.200 0.147 0.186 0.176 MSE( M e , t ) 0.063 0.067 0.097 0.063 0.068 0.119 0.073 0.092 0.078 0.098 0.082 Open in a new tab a The four paired one-sided Wilcoxon signed-rank tests on R 2 and MSE between M e,t and M e or M t all return p -value < 0.001, thus the null hypotheses are rejected, indicating model M e,t is significantly better than M e or M t . Limitation The manuscript focuses on providing algorithms for feature generation for proteins in general and is not restricted by protein size or type. However, we do notice some limitations for the proposed methods, as listed below. (1) The electrostatic features primarily focus on the protein electrostatic profile, therefore, they have limited impact on capturing the interplay between the protein and its surrounding environment. This limitation is reflected in the simulation results, where Coulombic energy prediction performs better than solvation energy prediction. In future work, we plan to enhance the electrostatic features by incorporating reaction potentials generated from the GB model at atomic centers to better characterize protein–solvent interactions. (2) When the dataset size is small, for the same types of features, the advantages of DNN over gradient boosting or random forest are difficult to justify. (3) Identifying the optimal values of p and L for generating electrostatic features is challenging, as they depend on factors such as dataset size, response types, and the performance of topological features. This limitation requires multiple tests to determine the optimal p and L during model training. We are planning to implement and test the Barycentric treecode approach as explained in the Supporting Information for potential solutions. Conclusion The project provides a novel approach for generating combined topological and electrostatic features to represent protein structures and charge distributions in a physically informed multiscale and uniform manner. Using these features, we can predict important protein properties, such as solvation energy and Coulomb energy. The numerical simulations validate our three main findings: larger dataset, higher-resolution electrostatic features, and the combination of features all improve the performance of the ML model. Beyond the prediction of solvation and Coulomb energies, the resulting electrostatic and topological features have the potential for broad applicability, particularly in representing proteins as uniform and multiscale features for ML models. In the future work, we will further investigate the role of the reaction potential, which will be computed using the much faster GB model, as discussed in the manuscript. We will also explore the use of the Barycentric Lagrange treecode for generating electrostatic features, in comparison with the Cartesian treecode. Supplementary Material ci5c03104_si_001.pdf (853.6KB, pdf) Acknowledgments This work of ES, MAT, and WG was supported in part by National Science Foundation (NSF) grants DMS-2110922 and DMS-2110869. ES was also partially supported by the NSF RTG-1840260 grant. We thank the SMU O’Donnell Data Science and Research Computing Institute (ODSRCI) for providing computing hardware and software. These resources together made this project possible. The code and data associated with this project are freely and publicly available in the GitHub repository: https://github.com/atalha-cmd/Biophysics_DNN , maintained by SMU graduate student Md Abu Talha. It contains the Fortran code for generating electrostatic features using Cartesian treecode, the Python code to generate topological features (a modified version as used in), the Python code for the neural network originating from Elyssa Sliheet’s thesis work, and the data required to reproduce the reported results. The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jcim.5c03104 . Additional information on models, methods, and results which are less essential but could assist the readers in better understanding or reproducing the proposed work ( PDF ) ‡. Vistra Corp The authors declare no competing financial interest. References Bragg W. L.. The Specular Reflection of X-rays. Nature. 1912;90:410. doi: 10.1038/090410b0. [ DOI ] [ Google Scholar ] Crowfoot D.. X-Ray Single Crystal Photographs of Insulin. Nature. 1935;135:591–592. doi: 10.1038/135591a0. [ DOI ] [ Google Scholar ] Aue W. P., Ernst R. R., Bartholdi E.. Two-dimensional spectroscopyApplication to nuclear magnetic resonance. J. Chem. Phys. 1976;64:2229–2246. doi: 10.1063/1.432450. [ DOI ] [ Google Scholar ] Cheng Y., Grigorieff N., Penczek P. A., Walz T.. A Primer to Single-Particle Cryo-Electron Microscopy. Cell. 2015;161:438–449. doi: 10.1016/j.cell.2015.03.050. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Berman H. M., Westbrook J., Feng Z., Gilliland G., Bhat T. N., Weissig H., Shindyalov I. N., Bourne P. E.. The Protein Data Bank. Nucleic Acids Res. 2000;28:235–242. doi: 10.1093/nar/28.1.235. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Benson D. A., Cavanaugh M., Clark K., Karsch-Mizrachi I., Lipman D. J., Ostell J., Sayers E. W.. GenBank. Nucleic Acids Res. 2012;41:D36–D42. doi: 10.1093/nar/gks1195. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Bepler T., Berger B.. Learning the protein language: Evolution, structure, and function. Cell Syst. 2021;12:654–669e3. doi: 10.1016/j.cels.2021.05.017. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Weissenow K., Rost B.. Are protein language models the new universal key? Curr. Opin. Struct. Biol. 2025;91:102997. doi: 10.1016/j.sbi.2025.102997. [ DOI ] [ PubMed ] [ Google Scholar ] Xu S., Onoda A.. Accurate and Rapid Prediction of Protein pKa: Protein Language Models Reveal the Sequence-pKa Relationship. J. Chem. Theory Comput. 2025;21:3752–3764. doi: 10.1021/acs.jctc.4c01288. [ DOI ] [ PubMed ] [ Google Scholar ] Lin Z., Akin H., Rao R.. et al. Evolutionary-scale Prediction of Atomic-level Protein Structure with a Language Model. Science. 2023;379:1123–1130. doi: 10.1126/science.ade2574. [ DOI ] [ PubMed ] [ Google Scholar ] Madani A., Krause B., Greene E. R.. et al. Large Language Models Generate Functional Protein Sequences Across Diverse Families. Nat. Biotechnol. 2023;41:1099–1106. doi: 10.1038/s41587-022-01618-2. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Zhang Z., Wayment-Steele H. K., Brixi G., Wang H., Kern D., Ovchinnikov S.. Protein Language Models Learn Evolutionary Statistics of Interacting Sequence Motifs. Proceedings of the National Academy of Sciences U.S.A. 2024;121:e2406285121. doi: 10.1073/pnas.2406285121. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Abramson J.. et al. Accurate structure prediction of biomolecular interactions with AlphaFold 3. Nature. 2024;630:493–500. doi: 10.1038/s41586-024-07487-w. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Gligorijević V., Renfrew P. D., Kosciolek T., Leman J. K., Berenberg D., Vatanen T., Chandler C., Taylor B. C., Fisk I. M., Vlamakis H., Xavier R. J., Knight R., Cho K., Bonneau R.. Structure-based protein function prediction using graph convolutional networks. Nat. Commun. 2021;12:3168. doi: 10.1038/s41467-021-23303-9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Sagendorf J. M., Mitra R., Huang J., Chen X. S., Rohs R.. Structure-based Prediction of Protein-nucleic Acid Binding Using Graph Neural Networks. Biophys. Rev. 2024;16:297–314. doi: 10.1007/s12551-024-01201-w. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Mitra R., Cohen A. S., Tang W. Y.. et al. RNAproDB: A Webserver and Interactive Database for Analyzing Protein–RNA Interactions. J. Mol. Biol. 2025;437:169012. doi: 10.1016/j.jmb.2025.169012. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Wang B., Zhao Z., Nguyen D. D., Wei G.-W.. Feature Functional Theory - Binding Predictor (FFT-BP) for the Blind Prediction of Binding Free Energies. Theor. Chem. Acc. 2017;136:55. doi: 10.1007/s00214-017-2083-1. [ DOI ] [ Google Scholar ] Wee J., Xia K.. Persistent Spectral Based Ensemble Learning (PerSpect-EL) for Protein-Protein Binding Affinity Prediction. Briefings Bioinf. 2022;23:bbac024. doi: 10.1093/bib/bbac024. [ DOI ] [ PubMed ] [ Google Scholar ] Townsend J., Micucci C. P., Hymel J. H., Maroulas V., Vogiatzis K. D.. Representation of molecular structures with persistent homology for machine learning applications in chemistry. Nat. Commun. 2020;11:3230. doi: 10.1038/s41467-020-17035-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Li P., Johnston H., Krasny R. A.. Cartesian Treecode for Screened Coulomb Interactions. J. Comput. Phys. 2009;228:3858–3868. doi: 10.1016/j.jcp.2009.02.022. [ DOI ] [ Google Scholar ] Cang Z., Wei G.-W.. TopologyNet: Topology based deep convolutional and multi-task neural networks for biomolecular property predictions. PLoS Comput. Biol. 2017;13:e1005690–27. doi: 10.1371/journal.pcbi.1005690. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Chen D., Chen Z., Chen C., Geng W., Wei G. W.. MIBPB: A software package for electrostatic analysis. J. Comput. Chem. 2011;32:756–770. doi: 10.1002/jcc.21646. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Achondo M. A., Chaudhry J. H., Cooper C. D.. An Investigation of Physics Informed Neural Networks to Solve the Poisson–Boltzmann Equation in Molecular Electrostatics. J. Chem. Theory Comput. 2025;21:3726–3744. doi: 10.1021/acs.jctc.4c01747. [ DOI ] [ PubMed ] [ Google Scholar ] Chen D., He C., Hu X., Mu L., Shen Z.. Deep Neural Network for Solving Poisson–Boltzmann Equations on Protein Surfaces. Journal of Machine Learning for Modeling and Computing. 2025;6:41–63. doi: 10.1615/JMachLearnModelComput.2025057332. [ DOI ] [ Google Scholar ] Huang Z., Li B., Wang Z., Zhang Z.. Neural network Poisson–Boltzmann electrostatics for biomolecular interactions. J. Comput. Phys. 2026;544:114446. doi: 10.1016/j.jcp.2025.114446. [ DOI ] [ Google Scholar ] Röcken S., Burnet A. F., Zavadlav J.. Predicting solvation free energies with an implicit solvent machine learning potential. J. Chem. Phys. 2024;161:234101. doi: 10.1063/5.0235189. [ DOI ] [ PubMed ] [ Google Scholar ] Ferraz-Caetano J., Teixeira F., Cordeiro M. N. D. S.. Explainable Supervised Machine Learning Model To Predict Solvation Gibbs Energy. J. Chem. Inf. Model. 2024;64:2250–2262. doi: 10.1021/acs.jcim.3c00544. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Case D. A.. et al. AmberTools. J. Chem. Inf. Model. 2023;63:6183–6191. doi: 10.1021/acs.jcim.3c01153. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Richards F. M.. Areas, Volumes, Packing, and Protein Structure. Annu. Rev. Biophys. Bioeng. 1977;6:151–176. doi: 10.1146/annurev.bb.06.060177.001055. [ DOI ] [ PubMed ] [ Google Scholar ] Connolly M. L.. Depth buffer algorithms for molecular modeling. J. Mol. Graphics. 1985;3:19–24. doi: 10.1016/0263-7855(85)80009-6. [ DOI ] [ Google Scholar ] Baker N. A.. Improving implicit solvent simulations: a Poisson-centric view. Curr. Opin. Struct. Biol. 2005;15:137–143. doi: 10.1016/j.sbi.2005.02.001. [ DOI ] [ PubMed ] [ Google Scholar ] Gilson M. K., Davis M. E., Luty B. A., McCammon J. A.. Computation of electrostatic forces on solvated molecules using the Poisson-Boltzmann equation. J. Phys. Chem. 1993;97:3591–3600. doi: 10.1021/j100116a025. [ DOI ] [ PubMed ] [ Google Scholar ] Geng W., Wei G. W.. Multiscale molecular dynamics using the matched interface and boundary method. J. Comput. Phys. 2011;230:435–457. doi: 10.1016/j.jcp.2010.09.031. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Chern I. L., Liu J.-G., Wang W. C.. Accurate Evaluation of Electrostatics for Macromolecules in Solution. Methods Appl. Anal. 2003;10:309–328. doi: 10.4310/maa.2003.v10.n2.a9. [ DOI ] [ Google Scholar ] LeVeque R., Li Z.. The Immersed Interface Method for Elliptic Equations with Discontinuous Coefficients and Singular Sources. SIAM J. Numer. Anal. 1994;31:1019–1044. doi: 10.1137/0731054. [ DOI ] [ Google Scholar ] Qiao, Z. ; Li, Z. ; Tang, T. . Finite difference scheme for solving the nonlinear Poisson-Boltzmann equation modeling charged spheres. J. Comput. Math. 2006, 24. [ Google Scholar ] Zhou Y. C., Zhao S., Feig M., Wei G. W.. High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources. J. Comput. Phys. 2006;213:1–30. doi: 10.1016/j.jcp.2005.07.022. [ DOI ] [ Google Scholar ] Yu S., Geng W., Wei G. W.. Treatment of geometric singularities in implicit solvent models. J. Chem. Phys. 2007;126:244108. doi: 10.1063/1.2743020. [ DOI ] [ PubMed ] [ Google Scholar ] Geng W., Yu S., Wei G. W.. Treatment of charge singularities in implicit solvent models. J. Chem. Phys. 2007;127:114106. doi: 10.1063/1.2768064. [ DOI ] [ PubMed ] [ Google Scholar ] Geng W., Zhao S.. A two-component Matched Interface and Boundary (MIB) regularization for charge singularity in implicit solvation. J. Comput. Phys. 2017;351:25–39. doi: 10.1016/j.jcp.2017.09.026. [ DOI ] [ Google Scholar ] Zauhar R. J., Morgan R. S.. The rigorous computation of the molecular electric potential. J. Comput. Chem. 1988;9:171–187. doi: 10.1002/jcc.540090209. [ DOI ] [ Google Scholar ] Juffer A., E B., van Keulen B., van der Ploeg A., Berendsen H.. The electric potential of a macromolecule in a solvent: a fundamental approach. J. Comput. Phys. 1991;97:144–171. doi: 10.1016/0021-9991(91)90043-k. [ DOI ] [ Google Scholar ] Liang J., Subranmaniam S.. Computation of molecular electrostatics with boundary element methods. Biophys. J. 1997;73:1830–1841. doi: 10.1016/s0006-3495(97)78213-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Boschitsch A. H., Fenley M. O., Zhou H.-X.. Fast Boundary Element Method for the Linear Poisson-Boltzmann Equation. J. Phys. Chem. B. 2002;106:2741–2754. doi: 10.1021/jp013607q. [ DOI ] [ Google Scholar ] Bordner A. J., Huber G. A.. Boundary element solution of the linear Poisson-Boltzmann equation and a multipole method for the rapid calculation of forces on macromolecules in solution. J. Comput. Chem. 2003;24:353–367. doi: 10.1002/jcc.10195. [ DOI ] [ PubMed ] [ Google Scholar ] Lu B., Cheng X., McCammon J. A.. A New-version-fast-multipole-method-accelerated electrostatic calculations in biomolecular systems. J. Comput. Phys. 2007;226:1348–1366. doi: 10.1016/j.jcp.2007.05.026. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Geng W., Krasny R.. A treecode-accelerated boundary integral Poisson-Boltzmann solver for electrostatics of solvated biomolecules. J. Comput. Phys. 2013;247:62–78. doi: 10.1016/j.jcp.2013.03.056. [ DOI ] [ Google Scholar ] Chen J., Tausch J., Geng W. A.. Cartesian FMM-accelerated Galerkin boundary integral Poisson-Boltzmann solver. J. Comput. Phys. 2023;478:111981. doi: 10.1016/j.jcp.2023.111981. [ DOI ] [ Google Scholar ] Barnes J., Hut P.. A hierarchical O­(NlogN) force-calculation algorithm. Nature. 1986;324:446–449. doi: 10.1038/324446a0. [ DOI ] [ Google Scholar ] Greengard L., Rokhlin V.. A fast algorithm for particle simulations. J. Comput. Phys. 1987;73:325–348. doi: 10.1016/0021-9991(87)90140-9. [ DOI ] [ Google Scholar ] Tausch J.. The Fast Multipole Method for arbitary Green’s Functions. Contemp. Math. 2003;329:307–314. doi: 10.1090/conm/329/05866. [ DOI ] [ Google Scholar ] Chen J., Geng W., Wei G.-W.. MLIMC: Machine learning-based implicit-solvent Monte Carlo. Chin. J. Chem. Phys. 2021;34:683–694. doi: 10.1063/1674-0068/cjcp2109150. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Chen J., Xu Y., Yang X., Cang Z., Geng W., Wei G.-W.. Poisson-Boltzmann-based machine learning model for electrostatic analysis. Biophys. J. 2024:2807–2814. doi: 10.1016/j.bpj.2024.02.008. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Edelsbrunner, H. ; Letscher, D. ; Zomorodian, A. . Topological Persistence and Simplification Proceedings 41st Annual Symposium on Foundations of Computer Science (FOCS); IEEE Computer Society, 2000.. [ Google Scholar ] Zomorodian A., Carlsson G.. Computing Persistent Homology. Discrete Comput. Geom. 2005;33:249–274. doi: 10.1007/s00454-004-1146-y. [ DOI ] [ Google Scholar ] Ghrist R.. Barcodes: The Persistent Topology of Data. Bull. Am. Math. Soc. 2008;45:61–75. doi: 10.1090/S0273-0979-07-01191-3. [ DOI ] [ Google Scholar ] Boissonnat J.-D., Maria C.. The Simplex Tree: An Efficient Data Structure for General Simplicial Complexes. Algorithmica. 2014;70:406–427. doi: 10.1007/s00453-014-9887-3. [ DOI ] [ Google Scholar ] Adams H., Emerson T., Kirby M., Neville R., Peterson C., Shipman P., Chepushtanova S., Hanson E., Motta F., Ziegelmeier L.. Persistence Images: A Stable Vector Representation of Persistent Homology. J. Mach. Learn. Res. 2017;18:1–35. [ Google Scholar ] Dolinsky T. J., Czodrowski P., Li H., Nielsen J. E., Jensen J. H., Klebe G., Baker N. A.. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. Nucleic Acids Res. 2007;35:W522–W525. doi: 10.1093/nar/gkm276. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] http://www.rcsb.org/pdb/home/home.do . Perlman D. A., Case D. A., Caldwell J. W., Ross W. S., Cheatham T. E., Debolt S., Ferguson D., Seibel G., Kollman P.. AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 1995;91:1–41. doi: 10.1016/0010-4655(95)00041-d. [ DOI ] [ Google Scholar ] Brooks B. R., Bruccoleri R. E., Olafson B. D., States D., Swaminathan S., Karplus M.. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983;4:187–217. doi: 10.1002/jcc.540040211. [ DOI ] [ Google Scholar ] Zhang C., Lu C., Jing Z., Wu C., Piquemal J.-P., Ponder J. W., Ren P.. AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids. J. Chem. Theory Comput. 2018;14:2084–2108. doi: 10.1021/acs.jctc.7b01169. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Su M., Yang Q., Du Y., Feng G., Liu Z., Li Y., Wang R.. Comparative Assessment of Scoring Functions: The CASF-2016 Update. J. Chem. Inf. Model. 2019;59:895–913. doi: 10.1021/acs.jcim.8b00545. [ DOI ] [ PubMed ] [ Google Scholar ] Pahari S., Sun L., Alexov E.. PKAD: a database of experimentally measured pKa values of ionizable groups in proteins. Database: Journal of Biological Databases and Curation. 2019;2019(article ID baz024):1–7. doi: 10.1093/database/baz024. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Liu Z., Li Y., Han L., Li J., Liu J., Zhao Z., Nie W., Liu Y., Wang R.. PDB–wide collection of binding data: current status of the PDBbind database. Bioinformatics. 2015;31:405–412. doi: 10.1093/bioinformatics/btu626. [ DOI ] [ PubMed ] [ Google Scholar ] UniProt Consortium. UniProt: the Universal Protein Knowledgebase in 2023. Nucleic Acids Res. 2016;51:D158–D169. doi: 10.1093/nar/gkw1099. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Steinegger M., Söding J.. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat. Biotechnol. 2017;35:1026–1028. doi: 10.1038/nbt.3988. [ DOI ] [ PubMed ] [ Google Scholar ] Wang L., Krasny R., Tlupova S.. A kernel-independent treecode based on barycentric Lagrange interpolation. Comput. Phys. Commun. 2025;28:1415–1436. doi: 10.4208/cicp.OA-2019-0177. [ DOI ] [ Google Scholar ] Sliheet, E.

Predicting Biomolecular Properties and Interactions Using Numerical, Statistical and Machine Learning Methods, Ph.D. thesis; Southern Methodist University, 2024.. [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials ci5c03104_si_001.pdf (853.6KB, pdf) Data Availability Statement The code and data associated with this project are freely and publicly available in the GitHub repository: https://github.com/atalha-cmd/Biophysics_DNN , maintained by SMU graduate student Md Abu Talha. It contains the Fortran code for generating electrostatic features using Cartesian treecode, the Python code to generate topological features (a modified version as used in), the Python code for the neural network originating from Elyssa Sliheet’s thesis work, and the data required to reproduce the reported results. Articles from Journal of Chemical Information and Modeling are provided here courtesy of American Chemical Society ACTIONS View on publisher site PDF (7.2 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top