A framework of finite element procedures for the analysis of proteins

Highlights

• A finite element framework for the analysis of proteins in solvent is newly developed.

• The effect of solvent damping is considered.

• The computational cost of this method is not dependent on the molecular size.

• Through numerical examples, the effectiveness of the method is demonstrated.

Abstract

Large-scale, functional collective motions of proteins and their supra-molecular assemblies occur in a physiological solvent environment at finite temperatures. The solution of these motions with standard molecular dynamics algorithms is computationally hardly possible when considering macromolecules. Much research has focused on alternative approaches that use coarse-graining to model proteins, but mostly in vacuum. In this paper, we incorporate realistically the physical effects of solvent damping into the finite element model of proteins. The proposed framework is based on Brownian dynamics and shown to be effective. An important advantage of the approach is that the computational cost is not dependent on the molecular size, which makes the long-time simulation of macromolecules possible. Using the proposed procedure, we demonstrate the analysis of a macromolecule in solvent—an analysis that has not been achieved before and could not be performed with a molecular dynamics algorithm.

Introduction

Protein motions such as conformational changes and folding/unfolding, generally occur in a physiological solvent, that is, a viscous environment within cells. Hence, to solve for the dynamical behavior of a protein, both the protein and the solvent should ideally be modeled simultaneously, as in all-atom, explicit-solvent molecular dynamics [1]. However, in practice, the time-integration of the full set of governing equations of motion in molecular dynamics is computationally hardly feasible, in particular when large length scale and long time scale motions need to be considered with the effects of the solvent.

Hence, to simulate protein motions, various coarse-grained modeling approaches have been developed. These models can describe approximately important protein motions that are hardly accessible using a molecular dynamics simulation. For example, protein folding and unfolding have been investigated, respectively, using the lattice models to coarse-grain the spatial discretization [2], [3], [4] and the steered molecular dynamics procedure to coarse-grain the time discretization [5]. Also, the elastic network model for coarse-grained normal mode analysis has been used to solve for the change of flexibility of proteins in large deformations [6], [7], [8], [9]. However, the effects of the solvent on the motion and the flexibility of proteins have been ignored in the elastic network model, and therefore, the predicted in-vacuum frequencies do not correspond to realistic time-scales and the physical normal modes of the protein [10], [11].

On the other hand, the Brownian dynamics formalisms include the solvent effects implicitly. The formalisms can be used to simulate biomolecular motions on a computer substantially faster than the molecular dynamics techniques and with finite element procedures [12] open an avenue to significantly advance the field. In 1978, Ermak and McCammon [13] proposed a generalized algorithm to simulate the Brownian dynamics of N particles, where hydrodynamic interactions were described by a 3N × 3N diffusion tensor. In the Ermak-McCammon procedure, the tensor needs to be Cholesky-decomposed [12] at each step to compute random displacements, resulting in a computation-time scaling of O(N³). Over the past four decades, researchers have developed several approaches to reduce this computational cost [14], [15], [16], [17] in order to make the long time-scale Brownian dynamics simulations of large biomolecules feasible. For example, a Chebyshev-polynomial approach was proposed by Fixman [16] for the approximation of the square-root of the diffusion tensor, which results in a computational cost that scales with O(N^2.25) [18]. Also, as another alternative to the direct Cholesky-decomposition of the diffusion tensor, Geyer and Winter [17] proposed the Truncated Expansion Ansatz, which scales with O(N²) by truncating the expansion of the hydrodynamic multi-particle correlations as two-body contributions at the second order. Recently, based on Krylov subspaces, Ando and coworkers [14] proposed a new approach to approximately compute the random Brownian displacements with a computation time scaling of O(N²). As an alternative to approximating the square-root of the diffusion tensor in order to speed up the Brownian dynamics simulations, the tensor may be also kept unchanged for several sequential time-steps [19], [20], [21], [22] or throughout the Brownian dynamics simulation as in our own recent work on DNA nanostructures [23].

A variety of Brownian dynamics packages are already available for simulating the protein dynamics from SDA [24] and Browndye [25], which use rigid-body models of proteins, to UHBD [26], BD_Box [27], and Brownmove [28], which use flexible models. Here, coarse-grained Brownian dynamics simulations have been used to analyze protein motions by employing bead models [29], [30], [31]. However, these models are complicated, and more importantly, they lead to bead overlapping [32], require volume and viscosity corrections [33], [34], and ignore the presence of protein atoms between bead pairs [35]. Additionally, although solvent friction takes place on the surface of proteins [36], the bead models used in the Brownian dynamics simulations assume that the frictional forces act at the centers of the α-carbon atoms (representative atoms in the protein) of amino acids which are the building blocks of proteins.

Here, we propose a novel framework of finite element procedures for the analysis of proteins. In this framework we model the protein and solvent environment more realistically with the frictional forces applied directly on the protein surface and without any overlapping and any correction for the volume and viscosity. The friction matrix due to the solvent damping is computed by embedding a protein in a Stokes fluid and establishing an influence matrix. Due to the specific physics, we do not solve a nonlinear fluid-structure interaction problem, like performed in many other fields, see for example [12], [37], [38]. The interaction matrix is obtained as usual in finite element analyses [12], [39], but of course with the specific conditions encountered in the case here considered, as detailed below. The computational cost to obtain the friction matrix for the Brownian dynamics simulation using ADINA version 9.3 (ADINA R&D, Inc, Watertown, MA, USA) [40] is quite reasonable.

In the following sections, we first discuss how the stiffness, mass, and friction matrices are obtained for the Brownian dynamics simulation using the finite element method, and show how to calculate the diffusion coefficients, which define the translational and rotational mobility of proteins in the solvent, from the friction matrix. Then, we give results obtained using the proposed method considering a simple case for which analytical solutions are available, and compare results for actual proteins with experimental data. Diffusion coefficients calculated for 10 proteins of various molecular weights, ranging from 7 kDa to 233 kDa (with 1 kDa = 1 kilodalton = $1.6605402\times {10}^{-21}$ g) are provided. We also give more detailed results for the proteins Taq polymerase and Lysozyme obtained using our Brownian dynamics, finite element simulation framework. These illustrate that the solvent-damping effects can significantly alter the normal modes of proteins. Finally, we show considering the protein GroEL that our proposed framework can be used efficiently to solve for the response of large proteins when a molecular dynamics solution is not feasible.