Applied Mathematics and Computational Science Program
For more information visit: https://cemse.kaust.edu.sa/amcs
Recent Submissions

DTiGEMS+: drug–target interaction prediction using graph embedding, graph mining, and similaritybased techniques.(Journal of Cheminformatics, Springer Science and Business Media LLC, 20200702) [Article]In silico prediction of drug–target interactions is a critical phase in the sustainable drug development process, especially when the research focus is to capitalize on the repositioning of existing drugs. However, developing such computational methods is not an easy task, but is much needed, as current methods that predict potential drug–target interactions suffer from high falsepositive rates. Here we introduce DTiGEMS+, a computational method that predicts Drug–Target interactions using Graph Embedding, graph Mining, and Similaritybased techniques. DTiGEMS+ combines similaritybased as well as featurebased approaches, and models the identification of novel drug–target interactions as a link prediction problem in a heterogeneous network. DTiGEMS+ constructs the heterogeneous network by augmenting the known drug–target interactions graph with two other complementary graphs namely: drug–drug similarity, target–target similarity. DTiGEMS+ combines different computational techniques to provide the final drug target prediction, these techniques include graph embeddings, graph mining, and machine learning. DTiGEMS+ integrates multiple drug–drug similarities and target–target similarities into the final heterogeneous graph construction after applying a similarity selection procedure as well as a similarity fusion algorithm. Using four benchmark datasets, we show DTiGEMS+ substantially improves prediction performance compared to other stateoftheart in silico methods developed to predict of drugtarget interactions by achieving the highest average AUPR across all datasets (0.92), which reduces the error rate by 33.3% relative to the secondbest performing model in the stateoftheart methods comparison.

Bayesian inference of spatially varying Manning’s n coefficients in an idealized coastal ocean model using a generalized KarhunenLoève expansion and polynomial chaos(Ocean Dynamics, Springer Science and Business Media LLC, 20200701) [Article]Bayesian inference with coordinate transformations and polynomial chaos for a Gaussian process with a parametrized prior covariance model was introduced in Sraj et al. (Comput Methods Appl Mech Eng 298:205–228, 2016a) to enable and infer uncertainties in a parameterized prior field. The feasibility of the method was successfully demonstrated on a simple transient diffusion equation. In this work, we adopt a similar approach to infer a spatially varying Manning’s n field in a coastal ocean model. The idea is to view the prior on the Manning’s n field as a stochastic Gaussian field, expressed through a covariance function with uncertain hyperparameters. A generalized KarhunenLoève (KL) expansion, which incorporates the construction of a reference basis of spatial modes and a coordinate transformation, is then applied to the prior field. To improve the computational efficiency of the method proposed in Sraj et al. (Comput Methods Appl Mech Eng 298:205–228, 2016a), we propose to use two polynomial chaos expansions to (i) approximate the coordinate transformation and (ii) build a cheap surrogate of the largescale advanced circulation (ADCIRC) numerical model. These two surrogates are used to accelerate the Bayesian inference process using a Markov chain Monte Carlo algorithm. Water elevation data are inverted within an observing system simulation experiment framework, based on a realistic ADCIRC model, to infer the KL coordinates and hyperparameters of a reference 2D Manning’s field. Our results demonstrate the efficiency of the proposed approach and suggest that including the hyperparameter uncertainties greatly enhances the inferred Manning’s n field, compared with using a covariance with fixed hyperparameters.

Solving Acoustic Boundary Integral Equations Using High Performance Tile LowRank LU Factorization.(High Performance Computing, Springer International Publishing, 20200618) [Book Chapter]We design and develop a new high performance implementation of a fast direct LUbased solver using lowrank approximations on massively parallel systems. The LU factorization is the most timeconsuming step in solving systems of linear equations in the context of analyzing acoustic scattering from large 3D objects. The matrix equation is obtained by discretizing the boundary integral of the exterior Helmholtz problem using a higherorder Nyström scheme. The main idea is to exploit the inherent data sparsity of the matrix operator by performing local tilecentric approximations while still capturing the most significant information. In particular, the proposed LUbased solver leverages the Tile LowRank (TLR) data compression format as implemented in the Hierarchical Computations on Manycore Architectures (HiCMA) library to decrease the complexity of “classical” dense direct solvers from cubic to quadratic order. We taskify the underlying boundary integral kernels to expose finegrained computations. We then employ the dynamic runtime system StarPU to orchestrate the scheduling of computational tasks on shared and distributedmemory systems. The resulting asynchronous execution permits to compensate for the load imbalance due to the heterogeneous ranks, while mitigating the overhead of data motion. We assess the robustness of our TLR LUbased solver and study the qualitative impact when using different numerical accuracies. The new TLR LU factorization outperforms the stateoftheart dense factorizations by up to an order of magnitude on various parallel systems, for analysis of scattering from largescale 3D synthetic and real geometries.

Hierarchical matrix approximations for spacefractional diffusion equations(Computer Methods in Applied Mechanics and Engineering, Elsevier BV, 20200611) [Article]Space fractional diffusion models generally lead to dense discrete matrix operators, which lead to substantial computational challenges when the system size becomes large. For a state of size N, full representation of a fractional diffusion matrix would require O(N2) memory storage requirement, with a similar estimate for matrix–vector products. In this work, we present H2 matrix representation and algorithms that are amenable to efficient implementation on GPUs, and that can reduce the cost of storing these operators to O(N) asymptotically. Matrix–vector multiplications can be performed in asymptotically linear time as well. Performance of the algorithms is assessed in light of 2D simulations of space fractional diffusion equation with constant diffusivity. Attention is focused on smooth particle approximation of the governing equations, which lead to discrete operators involving explicit radial kernels. The algorithms are first tested using the fundamental solution of the unforced space fractional diffusion equation in an unbounded domain, and then for the steady, forced, fractional diffusion equation in a bounded domain. Both matrixinverse and pseudotransient solution approaches are considered in the latter case. Our experiments show that the construction of the fractional diffusion matrix, the matrix–vector multiplication, and the generation of an approximate inverse preconditioner all perform very well on a single GPU on 2D problems with N in the range 105 – 106. In addition, the tests also showed that, for the entire range of parameters and fractional orders considered, results obtained using the H2 approximations were in close agreement with results obtained using dense operators, and exhibited the same spatial order of convergence. Overall, the present experiences showed that the H2 matrix framework promises to provide practical means to handle largescale space fractional diffusion models in several space dimensions, at a computational cost that is asymptotically similar to the cost of handling classical diffusion equations.

The large time profile for HamiltonJacobiBellman equations(arXiv, 20200608) [Preprint]Here, we study the largetime limit of viscosity solutions of the Cauchy problem for secondorder HamiltonJacobiBellman equations with convex Hamiltonians in the torus. This largetime limit solves the corresponding stationary problem, sometimes called the ergodic problem. This problem, however, has multiple viscosity solutions and, thus, a key question is which of these solutions is selected by the limit. Here, we provide a representation for the viscosity solution to the Cauchy problem in terms of generalized holonomic measures. Then, we use this representation to characterize the largetime limit in terms of the initial data and generalized Mather measures. In addition, we establish various results on generalized Mather measures and duality theorems that are of independent interest.

Accurately Solving Physical Systems with Graph Learning(arXiv, 20200606) [Preprint]Iterative solvers are widely used to accurately simulate physical systems. These solvers require initial guesses to generate a sequence of improving approximate solutions. In this contribution, we introduce a novel method to accelerate iterative solvers for physical systems with graph networks (GNs) by predicting the initial guesses to reduce the number of iterations. Unlike existing methods that aim to learn physical systems in an endtoend manner, our approach guarantees longterm stability and therefore leads to more accurate solutions. Furthermore, our method improves the run time performance of traditional iterative solvers. To explore our method we make use of positionbased dynamics (PBD) as a common solver for physical systems and evaluate it by simulating the dynamics of elastic rods. Our approach is able to generalize across different initial conditions, discretizations, and realistic material properties. Finally, we demonstrate that our method also performs well when taking discontinuous effects into account such as collisions between individual rods. A video showing dynamic results of our graph learning assisted simulations of elastic rods can be found on the project website available at http://computationalsciences.org/publications/shao2020physicalsystemsgraphlearning.html .

Approximation of an optimal control problem for the timefractional FokkerPlanck equation(arXiv, 20200605) [Preprint]In this paper, we study the numerical approximation of a system of PDEs with fractional time derivatives. This system is derived from an optimal control problem for a timefractional FokkerPlanck equation with time dependent drift by convex duality argument. The system is composed by a timefractional backward HamiltonJacobiBellman and a forward FokkerPlanck equation and can be used to describe the evolution of probability density of particles trapped in anomalous diffusion regimes. We approximate Caputo derivatives in the system by means of L1 schemes and the Hamiltonian by finite differences. The scheme for the FokkerPlanck equation is constructed such that the duality structure of the PDE system is preserved on the discrete level. We prove the well posedness of the scheme and the convergence to the solution of the continuous problem.

Computational Drugtarget Interaction Prediction based on Graph Embedding and Graph Mining(ACM, 20200522) [Conference Paper]Identification of interactions of drugs and proteins is an essential step in the early stages of drug discovery and in finding new drug uses. Traditional experimental identification and validation of these interactions are still timeconsuming, expensive, and do not have a high success rate. To improve this identification process, development of computational methods to predict and rank likely drugtarget interactions (DTI) with minimum error rate would be of great help. In this work, we propose a computational method for (DrugTarget interaction prediction using Graph Embedding and graph Mining), DTiGEM. DTiGEM models identify novel DTIs as a link prediction problem in a heterogeneous graph constructed by integrating three networks, namely: drugdrug similarity, targettarget similarity, and known DTIs. DTiGEM combines different techniques, including graph embeddings (e.g., node2vec), graph mining (e.g., path scores between drugs and targets), and machine learning (e.g., different classifiers). DTiGEM achieves improvement in the prediction performance compared to other stateoftheart methods for computational prediction of DTIs on four benchmark datasets in terms of area under precisionrecall curve (AUPR). Specifically, we demonstrate that based on the average AUPR score across all benchmark datasets, DTiGEM achieves the highest average AUPR value (0.831), thus reducing the prediction error by 22.4% relative to the secondbest performing method in the comparison.

Scattering cancellation technique for acoustic spinning objects(Physical Review B, American Physical Society (APS), 20200520) [Article]The scattering cancellation technique (SCT) has proved to be an effective way to render static objects invisible to electromagnetic and acoustic waves. However, rotating cylindrical or spherical objects possess additional peculiar scattering features that cannot be canceled by regular SCTbased cloaks. Here, a generalized SCT theory to cloak spinning objects, and hide them from static observers, based on rotating shells with different angular velocity is discussed. This concept is analytically and numerically demonstrated in the case of cylinders, showing that the generalized SCT operates efficiently in making rotating objects appear static to an external observer. Our proposal extends the realm of SCT and brings it one step closer to its practical realization that involves moving objects.

Asynchronous computations for solving the acoustic wave propagation equation(The International Journal of High Performance Computing Applications, SAGE Publications, 20200519) [Article]<jats:p> The aim of this study is to design and implement an asynchronous computational scheme for solving the acoustic wave propagation equation with absorbing boundary conditions (ABCs) in the context of seismic imaging applications. While the convolutional perfectly matched layer (CPML) is typically used for ABCs in the oil and gas industry, its formulation further stresses memory accesses and decreases the arithmetic intensity at the physical domain boundaries. The challenges with CPML are twofold: (1) the strong, inherent data dependencies imposed on the explicit timestepping scheme render asynchronous time integration cumbersome and (2) the idle time is further exacerbated by the load imbalance introduced among processing units. In fact, the CPML formulation of the ABCs requires expensive synchronization points, which may hinder the parallel performance of the overall asynchronous time integration. In particular, when deployed in conjunction with the multicoreoptimized wavefront diamond temporal blocking (MWDTB) approach for the inner domain points, it results in a major performance slow down. To relax CPML’s synchrony and mitigate the resulting load imbalance, we embed CPML’s calculation into MWDTB’s inner loop and carry on the time integration with finegrained computations in an asynchronous, holistic way. This comes at the price of storing transient results to alleviate dependencies from critical data hazards while maintaining the numerical accuracy of the original scheme. Performance and scalability results on various x86 architectures demonstrate the superiority of MWDTB with CPML support against the standard spatial blocking on various grid sizes. To our knowledge, this is the first practical study that highlights the consolidation of CPML ABCs with asynchronous temporal blocking stencil computations. </jats:p>

Automated methods for the comparison of natural languages(Computing and Visualization in Science, Springer Science and Business Media LLC, 20200518) [Article]Starting from the general question, if there is a connection between the mathematical capabilities of a student and his native language, we aim at comparing natural languages with mathematical language quantitatively. In [20] we set up an approach to compare language structures using Natural Language Processors (NLP). However, difficulties arose with the quality of the structural analysis of the NLP used just comparing simple sentences in different but closely related natural languages. We now present a comparison of different available NLPs and discuss the results. The comparison confirms the results from [20], showing that current NLPs are not capable of analysing even simple sentences such that resulting structures between different natural languages can be compared.

Estimating and forecasting COVID19 attack rates and mortality(Cold Spring Harbor Laboratory, 20200515) [Preprint]<jats:p>{We describe a model for estimating past and current infections as well as future deaths due to the ongoing COVID19 pandemic. The model does not use confirmed case numbers and is based instead on recorded numbers of deaths and on the age specific population distribution. A regularized deconvolution technique is used to infer past infections from recorded deaths. Forecasting is based on a compartmental SIRtype model, combined with a probability distribution for the time from infection to death. The effect of nonpharmaceutical interventions (NPIs) is modelled empirically, based on recent trends in the death rate. The model can also be used to study counterfactual scenarios based on hypothetical NPI policies.</jats:p>

Splice2Deep: An ensemble of deep convolutional neural networks for improved splice site prediction in genomic DNA(Gene: X, Elsevier BV, 20200513) [Article]Background: The accurate identification of the exon/intron boundaries is critical for the correct annotation of genes with multiple exons. Donor and acceptor splice sites (SS) demarcate these boundaries. Therefore, deriving accurate computational models to predict the SS are useful for functional annotation of genes and genomes, and for finding alternative SS associated with different diseases. Although various models have been proposed for the in silico prediction of SS, improving their accuracy is required for reliable annotation. Moreover, models are often derived and tested using the same genome, providing no evidence of broad application, i.e. to other poorly studied genomes. Results: With this in mind, we developed the Splice2Deep models for SS detection. Each model is an ensemble of deep convolutional neural networks. We evaluated the performance of the models based on the ability to detect SS in Homo sapiens, Oryza sativa japonica, Arabidopsis thaliana, Drosophila melanogaster, and Caenorhabditis elegans. Results demonstrate that the models efficiently detect SS in other organisms not considered during the training of the models. Compared to the stateoftheart tools, Splice2Deep models achieved significantly reduced average error rates of 41.97% and 28.51% for acceptor and donor SS, respectively. Moreover, the Splice2Deep crossorganism validation demonstrates that models correctly identify conserved genomic elements enabling annotation of SS in new genomes by choosing the taxonomically closest model. Conclusions: The results of our study demonstrated that Splice2Deep both achieved a considerably reduced error rate compared to other stateoftheart models and the ability to accurately recognize SS in other organisms for which the model was not trained, enabling annotation of poorly studied or newly sequenced genomes. Splice2Deep models are implemented in Python using Keras API; the models and the data are available at https://github.com/SomayahAlbaradei/Splice_Deep.git.

PositivityPreserving Adaptive RungeKutta Methods(arXiv, 20200513) [Preprint]Many important differential equations model quantities whose value must remain positive or stay in some bounded interval. These bounds may not be preserved when the model is solved numerically. We propose to ensure positivity or other bounds by applying RungeKutta integration in which the method weights are adapted in order to enforce the bounds. The weights are chosen at each step after calculating the stage derivatives, in a way that also preserves (when possible) the order of accuracy of the method. The choice of weights is given by the solution of a linear program. We investigate different approaches to choosing the weights by considering adding further constraints. We also provide some analysis of the properties of RungeKutta methods with perturbed weights. Numerical examples demonstrate the effectiveness of the approach, including application to both stiff and nonstiff problems.

Wave modeling of a reefsheltered coastal zone in the Red Sea(Ocean Engineering, Elsevier BV, 20200501) [Article]The coastal areas of the Red Sea are characterized by shallow banks of fringing and barrier reefs that provide protection against coastal hazards and erosion by dissipating wave energy. This study investigates the wave climate and extremes of a reefprotected coastal zone in the Red Sea using a highresolution coupled wave and circulation model, ADCIRC + SWAN, configured on an unstructured grid forced with the meteorological fields from highresolution regional atmospheric model. Our simulations suggest that the relatively narrow offshore reefs with steep forereef slopes dissipate 40–50% of the wave energy propagating towards the shoreline, and this is more pronounced during extremes. The impact of the coupling on determining the wave climate is negligible, but is significant for storms with ~10 cm higher significant wave height (Hs) during the observed period. The backreef wave climatology computed from 30year model simulations shows that the mean Hs distribution is uniform throughout the year, and extremes occur more often from February to May. Different return levels of Hs in the sheltered areas are estimated using extreme value analysis. Our results emphasize that preserving the complex offshore reefs is crucial for mitigating the coastal hazards of highenergy waves which are projected to increase with climate change.

3d Modeling and simulation of a harpsichord(Computing and Visualization in Science, Springer Science and Business Media LLC, 20200429) [Article]The mathematical characterization of the sound of a musical instrument still follows Schumann’s laws (Schumann in Physik der klangfarben, Leipzig, 1929). According to this theory, the resonances of the instrument body, “the formants”, filter the oscillations of the sound generator (e.g., strings) and produce the characteristic “timbre” of an instrument. This is a strong simplification of the actual situation. It applies to a point source and can be easily performed by a loudspeaker, disregarding the three dimensional structure of music instruments. To describe the effect of geometry and material of the instruments, we set up a 3d model and simulate it using the simulation system UG4 (Vogel et al. in Comput Vis Sci 16(4):165–179, 2013; Reiter et al. in Comput Vis Sci 16(4):151–164, 2014). We aim to capture the oscillation behavior of eigenfrequencies of a harpsichord soundboard and investigate how well a model for the oscillation behavior of the soundboard approximates the oscillation behavior of the whole instrument. We resolve the complicated geometry by several unstructured 3d grids and take into account the anisotropy of wood. The oscillation behavior of the soundboard is modeled following the laws of linear orthotropic elasticity with homogenous boundary conditions. The associated eigenproblem is discretized using FEM and solved with the iterative PINVIT method using an efficient GMG preconditioner (Neymeyr in A hierarchy of preconditioned eigensolvers for elliptic differential operators. Habilitation dissertation, University of Tübingen, 2001). The latter allows us to resolve the harpsichord with a high resolution hybrid grid, which is required to capture fine modes of the simulated eigenfrequencies. We computed the first 16 eigenmodes and eigenfrequencies with a resolution of 1.8 billion unknowns each on Shaheen II supercomputer (https://www.hpc.kaust.edu.sa/content/shaheenii). To verify our results, we compare them with measurement data obtained from an experimental modal analysis of a real reference harpsichord.

Dynamic programming bicriteria combinatorial optimization(Discrete Applied Mathematics, Elsevier BV, 20200429) [Article]We introduce a model for combinatorial optimization problems based on the notion of a circuit. This circuit builds a set of elements for optimization from oneelement sets attached to input nodes. It uses the operation of union of sets and functional operations on elements extended to sets of elements. We design a dynamic programming algorithm based on this circuit which constructs the set of Pareto optimal points for the problem of bicriteria optimization of elements described by the circuit relative to two cost functions. We tested this approach on nine known combinatorial optimization problems. The problems related to the matrix chain multiplication, optimal paths in directed graphs, and binary search trees are considered in detail.

Mesh generation for thin layered domains and its application to parallel multigrid simulation of groundwater flow(Computing and Visualization in Science, Springer Science and Business Media LLC, 20200420) [Article]The generation of detailed three dimensional meshes for the simulation of groundwater flow in thin layered domains is crucial to capture important properties of the underlying domains and to reach a satisfying accuracy. At the same time, this level of detail poses high demands both on suitable hardware and numerical solver efficiency. Parallel multigrid methods have been shown to exhibit near optimal weak scalability for massively parallel computations of density driven flow. A fully automated parameterized algorithm for prism based meshing of coarse grids from height data of individual layers is presented. Special structures like pinch outs of individual layers are preserved. The resulting grid is used as a starting point for parallel mesh and hierarchy creation through interweaved projected refinement and redistribution. Efficiency and applicability of the proposed approach are demonstrated for a parallel multigrid based simulation of a realistic sample problem.

Evaluating Cumulus Parameterization Schemes for the Simulation of Arabian Peninsula Winter Rainfall(Journal of Hydrometeorology, American Meteorological Society, 20200417) [Article]This study investigates the sensitivity of winter seasonal rainfall over the Arabian Peninsula (AP) to different convective physical parameterization schemes using a high resolution WRF model. Three different parameterization schemes: KainFritch (KF), BettsMillerJanjic (BMJ), and GrellFreitas (GF) are used in winter simulations from 2001 to 2016. Results from seasonal simulations suggest that simulated AP winter rainfall with KF is in best agreement with observed rainfall in terms of spatial distribution and intensity. Higher spatial correlation coefficients and less biases with observations are also obtained with KF. In addition, the regional moisture transport, cloud distribution, and cloud microphysical responses are better simulated by KF. The AP lowlevel circulation, characterized by the Arabian Anticyclone, is well captured by KF and BMJ, but its position is displaced in GF. KF is further more successful at simulating the moisture distribution in the lower atmosphere and atmospheric water plumes in the middle troposphere. The higher skill of rainfall simulation with the KF (and to some extent BMJ) is attributed to a better representation of the Arabian Anticyclone and subtropical westerly jet, which guides the upper tropospheric synoptic transients and moisture. In addition, the vertical profile of diabatic heating from KF is in better agreement with the observations. Discrepancies in representing the diabatic heating profile by BMJ and GF show discrepancies in instability and in turn precipitation biases. Our results indicate that the selection of subgrid convective parameterization in a highresolution atmospheric model over the AP is an important factor for accurate regional rainfall simulations.

Phase Consistent Ecological Domain Adaptation(arXiv, 20200410) [Preprint]We introduce two criteria to regularize the optimization involved in learning a classifier in a domain where no annotated data are available, leveraging annotated data in a different domain, a problem known as unsupervised domain adaptation. We focus on the task of semantic segmentation, where annotated synthetic data are aplenty, but annotating real data is laborious. The first criterion, inspired by visual psychophysics, is that the map between the two image domains be phasepreserving. This restricts the set of possible learned maps, while enabling enough flexibility to transfer semantic information. The second criterion aims to leverage ecological statistics, or regularities in the scene which are manifest in any image of it, regardless of the characteristics of the illuminant or the imaging sensor. It is implemented using a deep neural network that scores the likelihood of each possible segmentation given a single unannotated image. Incorporating these two priors in a standard domain adaptation framework improves performance across the board in the most common unsupervised domain adaptation benchmarks for semantic segmentation.