Refine
Year of publication
Language
- English (97)
Keywords
- optimal control (11)
- differential-algebraic equation (4)
- finite elements (4)
- metastability (4)
- trajectory storage (4)
- Lyapunov exponent (3)
- Sacker-Sell spectrum (3)
- exponential dichotomy (3)
- geodesic finite elements (3)
- rare events (3)
Project
In this review, we intend to clarify the underlying ideas and the relations
between various multigrid methods ranging from subset decomposition,
to projected subspace decomposition and truncated multigrid.
In addition, we present a novel globally convergent inexact active set method
which is closely related to truncated multigrid. The numerical properties
of algorithms are carefully assessed by means of a degenerate problem and
a problem with a complicated coincidence set.
We present a new inexact nonsmooth Newton method for the solution
of convex minimization problems with piecewise smooth, pointwise
nonlinearities. The algorithm consists of a nonlinear smoothing
step on the fine level and a linear coarse correction.
Suitable postprocessing guarantees global convergence
even in the case of a single multigrid step for each linear subproblem.
Numerical examples show that the overall efficiency
is comparable to multigrid for similar linear problems.
In this paper, we propose and investigate numerical methods based on QR factorization for computing all or some Lyapunov or Sacker-Sell spectral intervals for
linear differential-algebraic equations.
Furthermore, a perturbation and error analysis for these methods is presented. We
investigate how errors in the data and in the numerical integration affect the
accuracy of the approximate spectral intervals. Although we need to integrate
numerically some differential-algebraic systems on usually very long
time-intervals, under certain assumptions, it is shown that the error of the
computed spectral intervals can be controlled by the local error of numerical
integration and the error in solving the algebraic constraint.
Some numerical examples are presented to illustrate the theoretical results.
Motivated by the analysis of passive control systems, we undertake a detailed perturbation analysis of Hamiltonian matrices that have eigenvalues on the imaginary axis. We construct minimal Hamiltonian perturbations that move and coalesce eigenvalues of opposite sign characteristic to form multiple eigenvalues with mixed sign characteristics, which are then moved from the imaginary axis to specific locations in the complex plane by small Hamiltonian perturbations. We also present a numerical method to compute upper bounds for the minimal perturbations that move all eigenvalues of a given Hamiltonian matrix outside a vertical strip along the imaginary axis.
The PSurface Library
(2010)
We describe psurface, a C++ library that allows to store and access piecewise linear mappings between simplicial surfaces in $\R^2$ and $\R^3$. These mappings are stored in a graph data structure and can be constructed explicitly, by projection, or by surface simplification. Piecewise linear maps can be used, e.g., to construct boundary
approximations for finite element grids, and grid intersections for domain decomposition methods. In computer graphics the mappings allow to build level-of-detail representations as well as texture- and bump maps. We document the data structures and algorithms used and show how \psurface is used in the numerical analysis framework Dune
and the visualization software Amira.
This paper is devoted to the numerical approximation of Lyapunov and Sacker-Sell spectral intervals for linear differential-algebraic equations (DAEs). The spectral analysis for DAEs is improved and the concepts of leading directions and solution subspaces associated with spectral intervals are extended to DAEs. Numerical methods
based on smooth singular value decompositions are introduced for computing all or only some spectral intervals and their associated leading directions. The numerical algorithms as well as implementation issues are discussed in detail and numerical examples are presented to illustrate the theoretical results.
This paper presents three different adaptive algorithms for eigenvalue problems
associated with non-selfadjoint partial differential
operators. The basis for the developed algorithms is a homotopy method.
The homotopy method starts from a well-understood selfadjoint problem,
for which well-established adaptive methods are available.
Apart from the adaptive grid refinement, the progress of the homotopy as
well as the solution of the iterative method are adapted to balance the contributions
of the different error sources.
The first algorithm balances the homotopy, discretization and approximation errors with respect
to a fixed step-size $\tau$ in the homotopy.
The second algorithm combines the adaptive step-size control for the homotopy
with an adaptation in space that ensures an error below a fixed tolerance $\varepsilon$.
The third algorithm allows the complete adaptivity in space,
homotopy step-size as well as the iterative algebraic eigenvalue solver.
All three algorithms are compared in numerical examples.
In this work we propose a general framework for the structured perturbation
analysis of several classes of structured matrix polynomials in homogeneous
form, including complex symmetric, skew-symmetric, even and odd matrix polynomials. We introduce structured backward errors for approximate eigenvalues and eigenvectors and we construct minimal structured perturbations such that an approximate eigenpair is an exact eigenpair of an appropriately perturbed matrix polynomial. This work extends previous work for the non-homogeneous case (we include infinite eigenvalues) and we show that the structured backward errors improve the known unstructured backward errors.
The paper proposes goal-oriented error estimation and mesh refinement
for optimal control problems with elliptic PDE constraints using the value
of the reduced cost functional as quantity of interest. Error representation,
hierarchical error estimators, and greedy-style error indicators are derived and
compared to their counterparts when using the all-at-once cost functional as
quantity of interest. Finally, the efficiency of the error estimator and generated
meshes are demonstrated on numerical examples.
We investigate geodesic finite elements for functions with values in a space of zero curvature, like a torus or the M\"obius strip. Unlike in the general case, a closed-form expression for geodesic finite element functions is then available. This simplifies computations, and allows us to prove optimal estimates for the interpolation error in 1d and 2d. We also show the somewhat surprising result that the discretization by Kirchhoff transformation of the Richards equation proposed by Berninger et al. is a discretization by geodesic finite elements in the manifold $\mathbb{R}$ with a special metric.
We propose a generalization of the Structured Doubling Algorithm (SDA) to compute invariant subspaces
of structured matrix pencils
that arise in the context of solving linear quadratic optimal control problems.
The new algorithm is
designed to attain better accuracy when the classical Riccati equation approach for the solution of the optimal control problem is not well suited because
the stable and unstable invariant subspaces are not well separated (due to eigenvalues near or on the imaginary
axis) or in the case when the Riccati solution does not exist at all. We analyze the convergence
of the method and compare the new method with the classical SDA algorithm as well as some recent structured QR-methods.
This paper presents concepts and implementation of the finite element toolbox Kaskade 7, a flexible C++ code for solving elliptic and parabolic PDE systems. Issues such as problem formulation, assembly and adaptivity are discussed at the example of optimal control problems. Trajectory compression for parabolic optimization problems is considered as a case study.
Deflated and augmented Krylov subspace methods: Basic Facts and a Breakdown-free deflated MINRES
(2011)
In this paper we consider deflation and augmentation techniques for accelerating
the convergence of Krylov subspace methods for the solution of nonsingular linear
algebraic systems. The two techniques are conceptually different from
preconditioning. Deflation "removes" certain parts from the operator, while
augmentation adds a subspace to the Krylov subspace. Both approaches have been
used in a variety of methods and settings. For Krylov subspace methods that
satisfy a (Petrov-) Galerkin condition we show that augmentation can in general
be achieved implicitly by projecting the residuals appropriately and correcting
the approximate solutions in a final step. In this context, we analyze known
methods to deflate CG, GMRes and MinRes. Our analysis reveals that the recently
proposed RMinRes method can break down. We show how such breakdowns can be
avoided by choosing a special initial guess, and we derive a breakdown-free
deflated MinRes method. In numerical experiments we study the properties of
different variants of MinRes analyzed in this paper.
One of the most challenging problems in dynamic concurrent multiscale
simulations is the reflectionless transfer of physical quantities between the
different scales. In particular, when coupling molecular dynamics and finite
element discretizations in solid body mechanics, often spurious wave reflections
are introduced by the applied coupling technique. The reflected waves are
typically of high frequency and are arguably of little importance in the domain
where the finite element discretization drives the simulation.
In this work, we provide an analysis of this phenomenon.
Based on the gained
insight, we derive a new coupling approach, which neatly separates high and low
frequency waves. Whereas low frequency waves are permitted to
bridge the scales, high frequency waves can be removed by applying damping techniques without affecting the coupled share of the solution. As a consequence, our new method almost completely eliminates unphysical wave reflections and deals in a consistent way with waves of arbitrary frequencies. The separation of
wavelengths is achieved by employing a discrete $L^2$-projection, which acts as a
low pass filter. Our coupling constraints enforce matching in the range of this projection. With respect to the numerical realization this approach
has the advantage of a small number of constraints, which is computationally
efficient. Numerical results in one and two dimensions confirm our theoretical
findings and illustrate the performance of our new weak coupling approach.
We investigate optimal elliptic
regularity (within the scale of Sobolev spaces) of anisotropic
div--grad operators in three dimensions at a multi-material vertex on
the Neumann boundary part of a polyhedral spatial domain. The
gradient of a solution to the corresponding elliptic PDE (in a
neighbourhood of the vertex) is integrable to an index greater than
three.
We consider Large Deformation Diffeomorphic Metric Mapping of general $m$-currents. After stating an optimization algorithm in the function space of admissable morph generating velocity fields, two innovative aspects in this framework are presented and numerically investigated: First, we spatially discretize the velocity field with conforming adaptive finite elements and discuss advantages of this new approach. Second, we directly compute the temporal evolution of discrete $m$-current attributes.
We consider the mechanical coupling of a geometrically exact Cosserat rod to a linear elastic continuum. The coupling conditions are formulated in the nonlinear rod configuration space. We describe a Dirichlet--Neumann algorithm for the coupled system, and use it to simulate the static stresses in a human knee joint, where the Cosserat rods are models for the ligaments.
In optimal control problems with nonlinear time-dependent 3D PDEs, full 4D discretizations are usually prohibitive due to the storage requirement. For this reason gradient and Newton type methods working on the reduced functional are often employed. The computation of the reduced gradient requires one solve of the state equation forward in time, and one backward solve of the adjoint equation. The state enters into the adjoint equation, again requiring the storage of a full 4D data set. We propose a lossy compression algorithm using an inexact but cheap predictor for the state data, with additional entropy coding of prediction errors. As the data is used inside a discretized, iterative algorithm, lossy compression
maintaining a certain error bound turns out to be sufficient.
We introduce geodesic finite elements as a conforming way to discretize partial differential equations for functions $v : \Omega \to M$, where $\Omega$ is an open subset of $\R^d$ and $M$ is a Riemannian manifold. These geodesic finite elements naturally generalize standard first-order
finite elements for Euclidean spaces. They also generalize the geodesic finite elements proposed for $d=1$ by the author. Our formulation is equivariant under isometries of $M$, and hence preserves objectivity of continuous problem formulations. We concentrate on partial differential equations that can be formulated as minimization problems. Discretization leads to algebraic minimization problems on product manifolds $M^n$. These can be solved efficiently
using a Riemannian trust-region method. We propose a monotone multigrid method to solve the constrained inner problems with linear multigrid speed. As an example we numerically compute harmonic maps from a domain in $\R^3$
to $S^2$.
Logical modeling of biological regulatory networks gives rise to a representation of the system's dynamics as a so-called state transition graph. Analysis of such a graph in its entirety allows for a comprehensive understanding of the functionalities and behavior of the modeled system. However, the size of the vertex set of the graph is exponential in the number of the network components making analysis costly, motivating development of reduction methods. In this paper, we present results allowing for a complete description of an asynchronous state transition graph of a Thomas network solely based on the analysis of the subgraph induced by certain extremal states. Utilizing this notion, we compare the behavior of a simple multi-valued network and a corresponding Boolean network and analyze the conservation of dynamical properties between them. Understanding the relation between such coarser and finer models is a necessary step towards meaningful network reduction as well as model refinement methods.
Mathematical modeling often helps to provide a systems perspective on gene regulatory networks. In particular, qualitative approaches are useful when detailed kinetic information is lacking. Multiple methods have been developed that implement qualitative information in different ways, e.g., in purely discrete or hybrid discrete/continuous models. In this paper, we compare the discrete asynchronous logical modeling formalism for gene regulatory networks due to R. Thomas with piecewise affine differential equation models.
We provide a local characterization of the qualitative dynamics of a piecewise affine differential equation model using the discrete dynamics of a corresponding Thomas model. Based on this result, we investigate the consistency of higher-level dynamical properties such as attractor characteristics and reachability. We show that although the two approaches are based on equivalent information, the resulting qualitative dynamics are different. In particular, the dynamics of the piecewise affine differential equation model is not a simple refinement of the dynamics of the Thomas model.
The authors propose a recycling MINRES scheme for a solution of subsequent self-adjoint linear systems as appearing, for example, in the Newton process for solving nonlinear equations. Ritz vectors are automatically extracted from one MINRES run and then used for self-adjoint deflation in the next. The method is designed to work with a preconditioner and arbitrary inner products. Numerical experiments with nonlinear Schrödinger equations indicate a substantial decrease in computation time when recycling is used.
We consider a shape implant design problem that arises in the context of facial surgery. We introduce a reformulation as an optimal control problem, where the control acts as a boundary force. The state is modelled as a minimizer of a polyconvex hyperelastic energy functional. We show existence of optimal solutions and derive - on a formal level - first order optimality conditions. Finally, preliminary numerical results are presented.
The pole condition approach for deriving transparent boundary conditions is extended to the time-dependent, two-dimensional case. Non-physical modes of the solution are identified by the position of poles of the solution's spatial Laplace transform in the complex plane. By requiring the Laplace transform to be analytic on some
problem dependent complex half-plane, these modes can be
suppressed. The resulting algorithm computes a finite number of coefficients of a series expansion of the Laplace transform, thereby providing an approximation to the exact boundary condition. The resulting error decays super-algebraically with the number of coefficients, so relatively few additional degrees of freedom are
sufficient to reduce the error to the level of the discretization error in the interior of the computational domain. The approach shows good results for the Schroedinger and the drift-diffusion equation
but, in contrast to the one-dimensional case, exhibits instabilities for the wave and Klein-Gordon equation. Numerical examples are shown that demonstrate the good performance in the former and the instabilities in the latter case.
Flux variability analysis (FVA) is an important tool to further analyze the results obtained by flux balance analysis (FBA) on genome-scale metabolic networks. Standard FVA may predict unbounded fluxes through some reactions in the network even if the nutrient uptake rate is bounded. These fluxes violate the second law of thermodynamics. They may be eliminated by extending flux variability analysis with thermodynamic constraints.
We present a new algorithm for efficient flux variability (and flux balance) analysis with thermodynamic constraints, suitable for analyzing genome-scale metabolic networks. We first show that flux balance analysis with thermodynamic constraints is NP-hard. Then we derive a theoretical tractability result, which can be applied to metabolic networks in practice. We use this result to develop a new constraint programming algorithm Fast-tFVA for fast flux variability analysis with thermodynamic constraints (tFVA). Computational comparisons with previous methods demonstrate the efficiency of the new method. For tFVA, a speed-up of factor 30-300 is achieved.
In an analysis of genome-scale metabolic networks in the BioModels database, we found that in 485 out of 716 networks additional irreversible or fixed reactions could be detected.
This article describes Fortran 77 subroutines for computing eigenvalues and invariant subspaces
of Hamiltonian and skew-Hamiltonian matrices. The implemented algorithms are based on orthogonal
symplectic decompositions, implying numerical backward stability as well as symmetry
preservation for the computed eigenvalues. These algorithms are supplemented with balancing and
block algorithms, which can lead to considerable accuracy and performance improvements. As a
by-product, an efficient implementation for computing symplectic QR decompositions is provided.
We demonstrate the usefulness of the subroutines for several, practically relevant examples.
Stewart's recently introduced Krylov-Schur algorithm
is a modification of the implicitly restarted Arnoldi algorithm which
employs reordered Schur decompositions to perform restarts and de-
ations in a numerically reliable manner. This paper describes a variant
of the Krylov-Schur algorithm suitable for addressing eigenvalue
problems associated with products of large and sparse matrices. It
performs restarts and de
ations via reordered periodic Schur decompositions
and, by taking the product structure into account, it is
capable to achieve qualitatively better approximations to the eigenvalues
of small magnitude.
The role of larger bulges in the QR algorithm is controversial. Large bulges are infamous
for having a strong, negative influence on the convergence of the implicit shifted QR algorithm.
This paper provides a new explanation of this shift blurring effect by connecting the computation of
the first column of the shift polynomial to the notoriously ill-conditioned pole placement problem.
To avoid shift blurring, modern variants of the QR algorithm employ chains of tightly coupled tiny
bulges instead of one large bulge. It turns out that larger bulges still play a positive role in these
variants; a slight increase of the bulge sizes often results in considerable performance improvements.
We investigate the condition number for a complex eigenvalue of a real matrix under
real perturbations. Based on an explicit formula, it is shown that this number is never
smaller than 1/
p
2 times the corresponding condition number with respect to complex
perturbations. This result can be generalized to the condition number of an arbitrary
complex-valued function under real perturbations. This extends to related condition
numbers.
We give an algorithm to compute N steps of a convolution quadrature approximation
to a continuous temporal convolution using only O(N logN) multiplications and O(logN) active
memory. The method does not require evaluations of the convolution kernel, but instead O(logN)
evaluations of its Laplace transform, which is assumed sectorial. The algorithm can be used for the
stable numerical solution with quasi-optimal complexity of linear and nonlinear integral and integrodifferential
equations of convolution type. In a numerical example we apply it to solve a subdiffusion
equation with transparent boundary conditions.
In this paper we propose a new finite element realization of the Perfectly Matched
Layer method (PML-method). Our approach allows to deal with a wide class of
polygonal domains and with certain types of inhomogeneous exterior domains.
Among the covered inhomogeneities are open waveguide structures playing an essential
role in integrated optics. We give a detailed insight into implementation
aspects. Numerical examples show exponential convergence behavior to the exact
solution with the thickness of the PML sponge layer.
Lyapunov and exponential dichotomy spectral theory is extended
from ordinary differential equations (ODEs) to nonautonomous
differential-algebraic equations (DAEs). By using orthogonal
changes of variables, the original DAE system is transformed into
appropriate condensed forms, for which concepts such as Lyapunov
exponents, Bohl exponents, exponential dichotomy and spectral
intervals of various kinds can be analyzed via the resulting
underlying ODE. Some essential differences between the spectral
theory for ODEs and that for DAEs are pointed out. Numerical
methods for computing the spectral intervals associated with
Lyapunov and Sacker-Sell (exponential dichotomy) spectra are
derived by modifying and extending those methods proposed for ODEs. Perturbation theory and error analysis are discussed, as
well. Finally, some numerical examples are presented to illustrate
the theoretical results and the properties of the numerical
methods.
A Generic Grid Interface for Parallel and Adaptive Scientific Computing. Part I: Abstract Framework
(2007)
We give a mathematically rigorous definition of a grid for algorithms solving
partial differential equations. Unlike previous approaches, our grids have a
hierarchical structure. This makes them suitable for geometric multigrid
algorithms and hierarchical local grid refinement. The description is also
general enough to include geometrically nonconforming grids. The definitions
in this article serve as the basis for an implementation of an abstract grid
interface as C++ classes in the DUNE.
In a companion paper [Matheon-Preprint 403] we introduced an abstract definition of a parallel and adaptive hierarchical grid for scientific computing. Based on this
definition we derive an efficient interface specification as a set of C++ classes.
This interface separates the applications from the grid data structures.
Thus, user implementations become independent of the underlying grid
implementation. Modern C++ template techniques are used to provide an
interface implementation without big performance losses.
The implementation is realized as part of the
software environment DUNE.
Numerical tests demonstrate the flexibility and the efficiency of our approach.
In many applications such as data compression, imaging or
genomic data analysis,
it is important to approximate a given $m\times n$ matrix $A$
by a matrix $B$ of rank at most $k$ which is much smaller than $m$ and $n$.
The best rank $k$ approximation can be determined via
the singular value decomposition
which, however, has prohibitively
high computational complexity and storage requirements
for very large $m$ and $n$.
We present an optimal least squares algorithm for computing a rank $k$
approximation to an $m\times n$ matrix $A$ by reading
only a limited number of rows and columns of $A$.
The algorithm has complexity $\mathcal O(k^2\max(m,n))$ and
allows to iteratively improve given rank $k$
approximations by reading additional rows and
columns of $A$. We also show how this approach can be extended
to tensors and present numerical results.
In this paper, we discuss stability properties of positive descriptor systems in the continuous-time as well as in the discrete-time case. We present different characterisations of positivity and establish generalised stability criteria for the case of positive descriptor systems. We show that if the spectral projector onto the right finite deflating subspace of the matrix pair $(E,A)$ is non-negative, then all stability criteria for standard positive systems take a comparably simple form in the positive descriptor case. Furthermore, we provide sufficient conditions that guarantee entry-wise non-negativity along with positive semi-definiteness of solutions of generalised projected Lyapunov equations. As an application of the framework established throughout this paper, we exemplarily generalise two criteria for the stability of two switched standard positive systems under arbitrary switching to the descriptor case.
Adjoint Broyden a la GMRES
(2007)
It is shown here that a compact storage implementation of a quasi-Newton
method based on the adjoint Broyden update reduces in the affine
case exactly to the well established GMRES procedure. Generally,
storage and linear algebra effort per step are small multiples of $n\cdot k$,
where $n$ is the number of variables and $k$ the number of steps taken
in the current cycle. In the affine case the storage is exactly $(n+k)\cdot k$
and in the nonlinear case the same bound can be achieved if adjoints,
i.e. transposed Jacobian-vector products are available. A
transposed-free variant that relies exclusively on Jacobian-vector
products (or possibly their approximation by divided differences)
requires roughly twice the storage and turns out to be somewhat slower
in our numerical experiments reported at the end.
Structured eigenvalue conditioning and backward error of a class of polynomial eigenvalue problems
(2007)
Characterisations of simple eigenvalues of complex matrix polynomials with *-even/odd and *-palindromic/antipalindromic structures that have the same normwise condition number with respect to structure preserving and arbitrary perturbations are obtained. Here * denotes either the transpose T or the conjugate tranpose *. In the process we obtain formulae for the normwise structured condition number of simple eigenvalues of T-palindromic/antipalindromic and *-even/odd polynomials. Moreover, conditions under which the normwise structured backward error of approximate eigenvalues of such polynomials is equal to the unstructured error are also derived. These lead to complete characterisations of approximate eigenvalues that have the same structured and unstructured backward errors for the *-even/odd and T-even/odd polynomials.
The paper considers the time integration of frictionless dynamical contact problems between viscoelastic bodies in the frame of the Signorini condition. Among the numerical integrators, interest focuses on the contact-stabilized Newmark method recently suggested by Deuflhard et al., which is compared to the classical Newmark method and an improved energy dissipative version due to Kane et al. In the absence of contact, any such variant is equivalent to the Störmer-Verlet scheme, which is well-known to have consistency order 2. In the presence of contact, however, the classical approach to discretization errors would not show consistency at all because of the discontinuity at the contact. Surprisingly, the question of consistency in the constrained situation has not been solved yet. The present paper fills this gap by means of a novel proof technique using specific norms based on earlier perturbation results due to the authors. The corresponding estimation of the local discretization error requires the bounded total variation of the solution. The results have consequences for the construction of an adaptive timestep control, which will be worked out subsequently in a forthcoming paper.
We present an extension module for the Dune system. This module, called dune-subgrid, allows to mark elements of another Dune hierarchical grid. The set of marked elements can then be accessed as a Dune grid in its own right. dune-subgrid is free software and is available for download. We describe the functionality and use of dune-subgrid, comment on its implementation, and give two example applications.
First, we show how dune-subgrid can be used for micro-FE simulations of trabecular bone. Then we present an algorithm that allows to use exact residuals for the adaptive solution of the spatial problems of time-discretized evolution equations.
We present a hierarchical a~posteriori error analysis for the
minimum value of the energy functional in symmetric obstacle
problems. The main result is that the energy of the exact solution
is, up to data oscillation, equivalent to an appropriate
hierarchical estimator. The proof of the main result does not
invoke any saturation assumption. Moreover, we prove an a
posteriori error estimate indicating that the estimator from
\cite{RHWHoppe_RKornhuber_1994a} is asymptotically reliable and we
give sufficient conditions for the validity of a saturation
assumption. Finally, we corroborate and complement our theoretical
results with numerical experiments.
We introduce geodesic finite elements as a new way to discretize
the nonlinear configuration space of a geometrically exact Cosserat rod.
These geodesic finite elements naturally generalize standard one-dimensional
finite elements to spaces of functions with values in a Riemannian manifold.
For the special orthogonal group, our approach reproduces the
interpolation formulas of [Crisfield/Jelenic:1999].
Geodesic finite elements are
conforming and lead to objective and path-independent problem formulations.
We introduce geodesic finite elements for general Riemannian manifolds,
discuss the relationship between geodesic finite elements and
coefficient vectors, and estimate the interpolation error.
Then we use them to find static equilibria of hyperelastic Cosserat rods.
Using the Riemannian trust-region algorithm of [Absil/Mahony/Sepulchre:2008]
we show numerically that the discretization error depends optimally on
the mesh size.
The identification of metastable conformations of molecules plays an
important role in computational drug design. One main difficulty is the
fact that the underlying dynamic processes take place in high dimensional
spaces. Although the restriction of degrees of freedom to a few dihedral
angles significantly reduces the complexity of the problem, the existing
algorithms are time-consuming. They are mainly based on the approximation
of a transfer operator by an extensive sampling of states according
to the Boltzmann distribution and short-time Hamiltonian dynamics simulations.
We present a method which can identify metastable conformations
without sampling the complete distribution. Our algorithm is based
on local transition rates and uses only pointwise information about the
potential energy surface. In order to apply the cluster algorithm PCCA+,
we compute a few eigenvectors of the rate matrix by the Jacobi-Davidson
method. Interpolation techniques are applied to approximate the thermodynamical
weights of the clusters. The concluding example illustrates
our approach for epigallocatechine, a molecule which can be described by
seven dihedral angles.
In order to compute the thermodynamic weights of the different metastable conformations
of a molecule, we want to approximate the molecule’s Boltzmann distribution in a reasonable
time. This is an essential issue in computational drug design. The energy landscape of active
biomolecules is generally very rough with a lot of high barriers and low regions. Many of the
algorithms that perform such samplings (e.g. the hybrid Monte Carlo method) have difficulties
with such landscapes. They are trapped in low-energy regions for a very long time and cannot
overcome high barriers. Moving from one low-energy region to another is a very rare event. For
these reasons, the distribution of the generated sampling points converges very slowly against
the thermodynamically correct distribution of the molecule.
The idea of ConfJump is to use a priori knowledge of the localization of low-energy regions
to enhance the sampling with artificial jumps between these low-energy regions. The artificial
jumps are combined with the hybrid Monte Carlo method. This allows the computation of
some dynamical properties of the molecule. In ConfJump, the detailed balance condition is
satisfied and the mathematically correct molecular distribution is sampled.
Biochemical interactions are determined by the 3D-structure of the involved components –
thus the identification of conformations is a key for many applications in rational drug design.
ConFlow is a new multilevel approach to conformational analysis with main focus on
completeness in investigation of conformational space.
In contrast to known conformational analysis, the starting point for design is a space-based
description of conformational areas. A tight integration of sampling and analysis leads to an
identification of conformational areas simultaneously during sampling. An incremental
decomposition of high-dimensional conformational space is used to guide the analysis. A new
concept for the description of conformations and their path connected components based on
convex hulls and Hypercubes is developed. The first results of the ConFlow application
constitute a ‘proof of concept’ and are further more highly encouraging. In comparison to
conventional industrial applications, ConFlow achieves higher accuracy and a specified
degree of completeness with comparable effort.
In this paper, we investigate the interconversion processes of the major flame retardant - 1,2,5,6,9,10-hexabromocyclododecane (HBCD) - by the means of statistical thermodynamics based on classical force-fields. Three ideas will be presented. First, the application of classical hybrid Monte-Carlo simulations for quantum mechanical processes will be justified. Second, the problem of insufficient convergence properties of hybrid Monte-Carlo methods for the generation of low temperature canonical ensembles will be solved by an interpolation approach. Furthermore, it will be shown how free energy differences can be used for a rate matrix computation. The results of our numerical simulations will be compared to experimental results.
Keywords: Markov process, molecular dynamics, rate matrix
This article deals with an efficient sampling of the stationary distribution
of dynamical systems in the presence of metastabilities. For such
systems, standard sampling schemes suffer from trapping problems and
critical slowing down. Starting multiple trajectories in different regions of
the sampling space is a promising way out. The different samplings represent
the stationary distribution locally very well, but are still far away
from ergodicity or from the global stationary distribution. We will show
how these samplings can be joined together in order to get one global
sampling of the stationary distribution.
For the treatment of equilibrated molecular systems in a heat bath we
propose a transition state theory that is based on conformation dynamics.
In general, a set-based discretization of a Markov operator P does not
preserve the Markov property. In this article, we propose a discretization
method which is based on a Galerkin approach. This discretization
method preserves the Markov property of the operator and can be interpreted
as a decomposition of the state space into (fuzzy) sets. The
conformation-based transition state theory presented here can be seen as
a first step in conformation dynamics towards the computation of essential
dynamical properties of molecular systems without time-consuming
molecular dynamics simulations.