A DNN Biophysics Model with Topological and Electrostatic Features
Elyssa Sliheet Md Abu Talha Weihua Geng



smu]Elyssa Sliheet
esliheet@smu.edu

smu]Md Abu Talha
atalha@smu.edu

smu]Weihua Gengcor1
wgeng@smu.edu

[cor1]Corresponding author


[smu]Department of Mathematics, Southern Methodist University, Dallas, TX 75275 USA




§ ABSTRACT 
In this project, we provide a deep-learning neural network (DNN) based biophysics model to predict protein properties. 
The model uses multi-scale and uniform topological and electrostatic features generated with protein structural information and force field, which governs the molecular mechanics.
The topological features are generated using the element specified persistent homology (ESPH) while the electrostatic features are fast computed using a Cartesian treecode.
These features are uniform in number for proteins with various sizes thus the broadly available protein structure database can be used in training the network. 
These features are also multi-scale thus the resolution and computational cost can be balanced by the users. 
The machine learning simulation on over 4000 protein structures 
shows the efficiency and fidelity of these features in representing the protein structure and force field for the predication of their biophysical properties such as electrostatic solvation energy. 
Tests on topological or electrostatic features alone and the combination of both showed the optimal performance when both features are used.
This model shows its potential as a general tool in assisting biophysical properties and function prediction for the broad biomolecules using data from both theoretical computing and experiments.  



 


  






Electrostatics;
Poisson-Boltzmann;
Interface methods;
treecode;
multipole methods




A DNN Biophysics Model with  
Topological and Electrostatic Features

    [
    September 9, 2024
=====================================================================













Highlights

- Multiscale and Uniform Features

- High efficiency using the MtM Translation 

- Verified using GB and PB Data 

- DNN Topological Features + Electrostatic Features


Assignment

Eylssa & Weihua: Algorithm and Coding for the electrostatic features 

Elyssa: Verify the features using GB and PB Solvation Energy Data 

Elyssa: Repeat the topological DNN and add the electrostatic features for simulation




§ INTRODUCTION

One of the overarching themes of biology is that structure determines function.



§ MODELS AND ALGORITHMS

We introduce the theories and algorithms involved in this work. First in a comparison fashion, we introduce the two most popular implicit solvation models, the Poisson-Boltzmann model and the Generalized Born model, which are used to generate the core biological property of our concern: the electrostatic solvation energy, as the energy it takes bymoving the solutes into the solvent. Following that, we introduce the persistent homology, which produce the topological features. Finally
we introduce the Catersian treecode, whose modification evolves into the algorithms for us to generate the electrostatic features.  




















 §.§ The Poisson-Boltzmann model
























Figure 1 depicts the popular implicit solvent models. 
In Fig. 1(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 Fig. 1(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 salt ions  <cit.>.
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 weak electrostatic potential
one obtains  the linearized Poisson-Boltzmann (PB) model <cit.>,


    _1^2 _1(x)    = 
    -∑_k=1^N_c q_k (x-x_k),    x∈_1, 
    _2^2 _2(x)    =
    ^2_2(x),    x∈_2, 
    _1(x)    = _2(x),_1
    _1 n( x) = _2_2 n( x), 
       x∈Γ,


where  is the inverse Debye length measuring the salt concentration,
and
the potential satisfies a zero far-field boundary condition.

The PB model governs the electrostatic potential ϕ in the entire space. Theoretically after ϕ
is obtained, its gradient will produce 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 Γ <cit.>). Our attention for this project is on the energy side as described below. 
The electrostatic potential energy is given as 

    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(ϕ_reac( x_k)+ϕ_coul( x_k)) = E_solv+E_coul


where ρ( x) = ∑_k=1^N_c q_k (x-x_k) is the charge density as a sum of partial charges weighted delta function and  the E_solv =  1/2∑_k=1^N_c q_k ϕ_reac( x_k) term is the solvation energy, the energy it takes for the protein to solvate from the vacuum to the solvent. The ϕ_reac is the reaction potential as the remaining component when Coulomb potential ϕ_coul is taken away from the total electrostatic potential ϕ. 


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,
(4) the domain is unbounded.















To overcome these numerical difficulties, 
several  finite difference interface methods have been developed.
By assuming that the interface is aligned with a mesh line,
a jump condition capture scheme has been developed in <cit.>.
Based on a Cartesian grid, the immersed interface method (IIM) <cit.>
has been applied to solve the PB equation 
<cit.>, in which the jump conditions can be rigorously enforced based on Taylor expansions. 
For the purpose of dealing with arbitrarily shaped  dielectric interfaces based
on a simple Cartesian grid,
a matched interface and boundary (MIB) PB solver
<cit.> 
has been developed through rigorous treatments of geometrical and charge singularities.





Boundary element methods (BEM) for the PB model were developed
later <cit.> 
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.
In the original BEMs, these advantages were offset by
the high cost of evaluating the interactions among the elements,
but fast summation schemes have been developed to reduce the 
cost <cit.>
and
in our previous work we employed a O(N log N) Cartesian treecode  <cit.> and later a O(N) Cartesian Fast Multipole Method (FMM) <cit.> for this purpose. 
In developing these boundary integral PB solvers, we investigate and resolved many interesting numerical challenges e.g. preconditioning of matrix whose conditioner number increases when triangulation quality is reduced, the parallelization of the treecode algorithm using MPI <cit.> and the parallelization of the boundary integral PB solver using GPU <cit.>. The developed solvers can effeciently help us to produce electrostatic potentials, which can be further used to compute protein properties such as binding energy <cit.> and pKa values <cit.>. 






 §.§ Generalized Born Model

The solvation energy can also be efficiently computed by using the Generalized Born's method, which is an approximation to the Poisson's model without considering the ions. 
This model serves to reduce the computational complexity of the PB model, which determining the solution requires solving a 3D partial differential equation.  





To derive the GB model, we first 
we adopt some classic results for electrostatics from Jackson's textbook <cit.>. Using the analogue from discrete point charge distribution to continuous charge density distribution,  the potential energy of a charged system with density distribution ρ( r) takes the form

    W = 1/2∫ρ( r) ϕ(  r) d r = 1/2∫ E( r) · D( r) d r

where the electric field E = -∇ϕ, and electric displacement D=ϵ E. 

Now consider assembling a charge to the center at the origin of a sphere with radius r_i. The sphere separates the domain with ϵ_in for r <a_i and ϵ_out for r <a_i. The assembly takes the energy 

    G_i = 1/8π∫ D· D/ϵ d r≈1/8π∫_r<a_iq_i^2/r^4ϵ_in dx + 1/8π∫_r>a_iq_i^2/r^4 ϵ_out dx

where the Coulomb field approximation D = q  r/r^3 is used. 

The solvation energy is the energy difference when ϵ for r>a_i changes from ϵ_in (unsolvated) to ϵ_out (solvated) given as

    Δ G_solv,i = 1/8π(1/ϵ_out - 1/ϵ_in ) ∫_r>a_iq_i^2/r^4dx
    
    
    = (1/ϵ_out - 1/ϵ_in) q_i^2/2a_i.

This result is consist with the Poisson-Boltzmann model of a solvated spherical cavity with centered charges as summarized in <cit.>, which is a special case from Kirkwood's derivation of a series of spherical harmonics for a spherical cavity containing arbitrary multiple charges <cit.>. 
If we treat the molecule as a collection of spherical atoms, using Eq. (<ref>) the total solvation energy is

    Δ G_solv =1/2(1/ϵ_out-1/ϵ_in)( ∑_i=1^N q_i^2/a_i + ∑_j ≠ i^N q_i q_j/r_ij)
    ≈1/2 (1/ϵ_out-1/ϵ_in) ∑_i=1^N ∑_j=1^Nq_i q_j/f_ij^GB

where the first equation has two terms: a sum of individual Born terms and pairwise Coulombic terms. 
Note the actual molecule is better represented by a dielectric interface 
(e.g. solvent excluded surface (SES)) which separates inside domain Ω_in and outside domain Ω_out. The charge q_i is located at the center r_i of a sphere with radius a_i.  We assume Ω_in contains all these spheres. Thus the model using dielectric interface requires a step further as been approximated in the second equation in Eq. (<ref>). 

Here f_ij^GB is the effective Born radii (i=j) or effective interaction distance (i  j). 
Assuming the Born Radii R_i's are obtained, a popular estimation of 
f_ij^GB is:

    f_ij^GB(r_ij) = (r_ij^2 +R_iR_j e^-r_ij^2/4R_iR_j)^1/2

where r_ij is the distance between the atomic centers of atoms i and j.
Note R_i depends not only on a_i, but also on radii and relative positions of all other atoms. 
The estimation of effective Born radii is an active research area. From Poisson-Boltzmann theory, the perfect Born radii is given as

    R_i = 1/2 (1/ϵ_out-1/ϵ_in) q_i/Δ G_PB,i

with Δ G_PB,i as the solvation energy from the PB model with all except ith charge muted. 
There are many approaches to approximate R_i and a typical one is the volume integration via quadrature as

    R_i^-1 = 1/4π∫_Ω_outq_i^2/r^4 d r = q_i^2/a_i-1/4π∫_Ω_in/B( r_i, a_i)q_i^2/r^4 d r.

Computing R_i using FFT leads to the computational cost of sovaltion energy to O(n^3log n+N+N^2) <cit.> against the cost of solving the PB model of O(n^6+Nn^2+N) with finite difference method (iterative solver, matrix structure not considered), where n is the number of grid point in x,y,and z direction, assuming cube-alike computational domain.






As shown in Fig. 1(c), we define the perfect Born Radius of the atom centered at x_k as the radius of the red sphere that makes the solvation energy of the red sphere 
is equal to the solvation energy of the protein 
when there is only one charge q_k located at x_k for both the sphere and the protein. 
The GB model calculates the solvation energy as a sum of interactions between spherical atoms of 
Born radii R_k for k=1,⋯,N_c

    E_solv = 1/2(1/ϵ_1-1/ϵ_2)∑_i,j=1^N_cq_i q_j/f_ij

with effective interaction distance
f_ij = √(r_ij^2+R_iR_jexp(-r^2_ij/4R_iR_j)).

The approximation of R_k efficiently and accurately is a hot research topic, aiming to obtain the perfect Born radii rapidly. We briefly mention the GB model here due to its efficiency to compute solvation energy compared with the cost to solve the PB model. Especially in the machine learning model, using quantities produced from GB can be a great choice of features when PB quantities need to be fast predicted <cit.>.  









  












































































































































































































































































































 §.§ Topological Features

The fundamental task of topological data analysis is to extract topological invariants as the 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 by algebraic or geometric way. 
In addition, topological invariants in a discrete data set 
can be studied using simplicial homology which uses a specific rule to identify simplicial complexes from simplexes. 
Here the simplex represents the simplest possible polytope 
in any given dimension like point, 
line segment, 
triangle, 
tetrahedron, etc. 
Furthermore, filtration and persistent homology 
can identify and connect complexes at different level of complexity 
and the appearance and disappearance of the homology group. 
This brings us an analog to the four levels of structure of proteins 
(Primary, Secondary, Tertiary, and Quaternary) 
and the change in protein structure during its folding pathway, as well as
binding affinity at different location and orientation. 

To make the topological features uniform and physically informed, 
we choose to use Element Specific Persistent Homology (ESPH) 
to extract topological features at different level of complexity. 
An example of Element Specific is that we can use the collection of vertices V_𝒫_k for the unit 𝒫_k in 𝒫={ CC, CN,CO,CS,CH, NN,⋯,SS,SH,HH}  as mentioned earlier to generate topological features.  
ESPH has been used for successfully predicting protein-ligand binding <cit.>. Once the set is determined and numerically calculated, the barcode can be generated using software GUDHI <cit.>, followed by the algorithm below to generate the vector of topological features. 

We define the collection of barcodes as 𝔹(α, 𝒞, 𝒟), where

α : atom labels (i.e., protein, ligand, mutated residue)

𝒞 : type of simplicial complex (i.e., Ribs or Cech)

𝒟 : dimension (i.e., Betti-0, Betti-1, 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 are defined as:


    V^b    = || {(b_j, d_j) ∈𝔹(α, 𝒞, 𝒟) | (i-1)L/n ≤ b_j ≤ iL/n } ||, 1≤ i < n, 
     V^d    = || {(b_j, d_j) ∈𝔹(α, 𝒞, 𝒟) | (i-1)L/n ≤ d_j ≤ iL/n } ||, 1≤ i < n, 
     V^p    = || {(b_j, d_j) ∈𝔹(α, 𝒞, 𝒟) | (i-1)L/n ≥ b_j, iL/n ≤ d_j } ||, 1≤ i < n,


These vectors can be viewed as (1D) images. Each pixel is associated with, m channels that describe different element type, mutation status, topological dimension, and topological event (birth and death).   




 §.§ Electrostatic Features 



  §.§.§ General introduction

In majority of the molecular simulations, including Monte Carlo simulation,  Brownian Dynamics,  Molecular Dynamics, etc, the electrostatic interactions are characterized by the interactions between the partial charges assigned at the atomic centers by the force field, which are determined by experiment or quantum chemistry.  




































The algorithm for obtaining the mulit-scale, physics-informed, uniform electrostatic features is explained in Fig. 
<ref>. In short, the electrostatic of the protein from q and ϕ_reac at the charge locations will be represented by the point-multipoles, whose moments will be calculated using treecode or FMM combined. To understand the point-multipoles, we provide some explanations below. 


Consider a protein with N_c atom and its multi-scale N_d point multipole representation as shown in Fig. <ref> for our machine learning model. For n=1,⋯, N_c, at the center of the nth atom, i.e., r_n=(x_n,y_n,z_n),  
 the nth permanent order 2 (for example) multipole M^n consists of 13 components:  
M^n = [ q^n, d^n_x,  d^n_y,  d^n_z, Q^n_xx, Q^n_xy, …, Q^n_zz]^T, where q, d_i, Q_ij for i,j=1,2,3 are the moments of the monopole, dipole, quadruple in suffix notation.
Using this notation, the permanent charge at  r_n can be written as  <cit.>

    ρ^n( r)=q^nδ( r- r_n)+d_i^n∂_iδ( r- r_n)
    	+Q^n_ij∂_ijδ( r- r_n),

A key idea is that the Coulomb potential G^n

governed by the Gauss's law - Δ G^n = 4 πρ^n in the free space is expressed in terms of the Green's function

    G^n( r)=1/| r- r_n|q^n + r_i - r_n,i/| r- r_n|^3d_i^n
    + (r_i - r_n,i)(r_j - r_n,j)/2| r- r_n|^5 Q^n_ij.

For all permanent multipoles M = [   M^1,  M^2, …,  M^N_d]^T, the total Coulomb potential is additive such as 
G^ M( r) = ∑_n=1^N_d G^n( r) by the superposition principle. 


 





For the point-multipole approach in the left picture of Fig <ref>, our goal is thus the computation of M accurately and efficiently. 
In fact, the computational cost is O(N_c) using the strategies in <cit.>, i.e. 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. These moments are intrinsic properties of the cluster thus can serve as features for the protein, which carries simplified and important information.  The number of features up to level L cluster is 

    N_f(p,L) = N_p(1+8+8^2+⋯ + 8^L) = N_p8^L+1-1/7

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 pth order multipole is used. 

Alternatively, for the Barycentric treecode approach in the right picture of Fig <ref>, the computational cost is O(N_clog N_c), 
and the number of terms in the similar sense of N_p is c_p^3 where c_p is the number of Chebyshev nodes in one direction <cit.>.  




  §.§.§ Implementation details


1. 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.

2. Two parameters p and L are given, where p is the order in the Cartesian Taylor expansion and L are 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) = N_p(1+8+8^2+⋯ + 8^L) = N_p8^L+1-1/7=p(p+1)(p+2)(8^L+1-1)/42.

These N_f numbers are ordered by level from 0 to L and the coordinates of cluster centers in each level (more details to be given). This determines the dimension of our final feature vector F ∈ℝ^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   18724 

2    10   90   730   5850    46810 

3    20   180   1460   11700    93620 

4    35   315   2555   20475   163835 




Steps:

1. Read in the particle location and charges

2. From order p and level L to calculate number of features N_f

3. Order the particles 

4. Build the tree with L levels; Note each cluster has variables of:






 §.§ The machine learning models


The proposed research has the following components. 

(a) A DNN based machine learning model using physics informed and multi-aspect features. These features are in  algebra,  topology, geometry, and electrostatics 
and are extracted using mathematical algorithms and their related software. 
These features are multi-scale and uniform. 
Being multi-scale, we mean depending on the accuracy requirement 
and available computing powers, 
the number of features can be adjusted as needed.
Being uniform, we mean once the scale is set, 
the number of features is fixed for proteins of different sizes. 

(b) Algorithms to produce electrostatic features using 
the partial charge distribution that is assigned by atomic location from X-ray or NMR and the force field.
We will consider both the pairwise Coulomb interaction between partial charges 
and the reaction potentials output 
from solving the Poisson-Boltzmann model. 

(c) Biological applications of designed machine learning model. 
We start with applications we had experience such as  electrostatic solvation energy, 
free energy and binding energy if protein and its ligand are both involved; 
protein pKa prediction; 
protein Monte Carlo simulation; and eventually protein molecular dynamics.   












 §.§ A DNN based machine learning model with algebraic, topological, geometric, and electrostatic features





The main goal of this project is to discover the hidden and useful information embedded 
in the protein structural data from the protein data bank in a simple and abstract way. 
The protein structure data and force field 
include bonded and non-bonded interactions with which 
the molecular simulation can be performed. 
We categorize these information into algebraic, topological, geometric, and electrostatic features. 
These four aspects are all important to the machine learning model. 
Here we lightly touch the first three  
and put our focus on the electrostatic features as our main research focus. 
Based on our two-decade experience in developing numerical algorithms in electrostatic interations, 
our proposed algorithms are novel and practical with the combination of both accuracy and efficiency.

The DNN based Machine Learning model is shown in Fig. <ref>, which uses the available protein structural data <cit.>, force field such as AMBER <cit.>, CHARMM<cit.>, AMOEBA <cit.>, etc. and known protein properties repositories  such as PDBbind Database <cit.>, Protein pKa Database <cit.>, etc to mathematically generate algebraic, topological, geometric, and electrostatic features to train a learned model 
and then use the learned model to predict unknown properties for protein with available structures. 






































































§ RESULTS



 §.§ Electrostatic Features Only



 §.§ Topological Features Only



 §.§ Electrostatic and Topological Features Combined




§ DATA AND SOFTWARE DISSEMINATION



Questions:

1. What's the purpose of Run_alpha_hydro()?

2. The reason that only 1000+ proteins can have both topological features and electrostatic features. 

The 4000+ proteins with both PB and GB solvation energy computed has only protein.pqr files. 

Solution: modify the files which generate topological features to only use pro.pqr file for atoms [C,N,O,S] 

3. For all the performance parameters: MSE, PCC, R^2, MAPE, and the scatter plot, are these results from testing data?

4. How much difference will the simulations at different time have?
 



In the present work, 
the selected 4294? protein structures are obtained from the PDBbind v2015 refined set and core set, and PDBbind v2018 refined set as the training set <cit.>. The collection has proteins sized  from 997 to 27,713 atoms.  


A data pre-processing is required before a PB solver can be called. 

The protein structures in the original data set are protein-ligand complexes. 
Missing atoms and side-chains are filled using the protein preparation wizard utility of the Schrodinger 2015-2 Suite with default parameter setting. 
The Amber ff14SB general force field is applied for the atomic van der Waals radii and partial charges. 



 §.§ Topology Features Generation

The following software needs to be installed.

https://gudhi.inria.fr/index.htmlGUDHI 


Currently a protein .pdb file and a ligand .mol2 are required as the inputs for the generation of topological features. 

1. Get_structure(pdb_file, mol2_file ,'complex.npz')

Finding atoms [C,N,O,S] in protein within cut-off distance 50Å from atoms [C,N,O,S,P,F,Cl,Br,I] in ligand, and save all these atoms in terms of position and types in the complex.npz file. 


2. Run_alpha('complex.npz', 'protein.PH', 'ligand.PH', 'complex.PH')

Use the atoms positions in protein, ligand, and complex to create Cech/Alpha complex in terms of barcode of simplexes and then convert the barcode information to Betti numbers.  


3. Run_alpha_hydro('complex.npz', 'protein_C-C.PH', 'ligand_C-C.PH', 'complex_C-C.PH')



4. PrepareData('complex.npz', 'protein.PH', 'complex.PH', 'protein_C-C.PH', 'complex_C-C.PH', 'complex_digit.npy')

Turn Betti numbers into Topological Features in terms of 1-D like images. 









































































































ieeetr




