Mathematical modeling in biology: Systems capable to live and die

The tendencies in the mathematical modeling of biological systems which under certain conditions are capable of supporting survival were reviewed. Its distinctive features were considered and inherent drawbacks analyzed. Models of limited life span in the animal organism were discussed, and a general scheme of the vital processes in organism, ageing, and death was presented. The oxidative theory-based homeostatic approach to ageing was shown to enable one to combine and generalize the existing theories of ageing. Consideration was given to the biotechnical applications of the homeostatic approach.     

Exact enumeration of three-dimensional lattice proteins

We present an algorithm for the exhaustive enumeration of all monomer sequences and conformations of short lattice proteins as described by the hydrophobic-polar (HP) model. The algorithm is used for an exact identification of all designing sequences of HP proteins consisting of up to 19 monomers whose conformations are represented by interacting self-avoiding walks on the simple cubic lattice. Employing a parallelized implementation on a Linux cluster, we generate the complete set of contact maps of such walks.     

An enhanced version of SMMP—open-source software package for simulation of proteins

We describe a revised and updated version of the program package SMMP (Simple Molecular Mechanics for Proteins) [F. Eisenmenger, U.H.E. Hansmann, Sh. Hayryan, C.-K. Hu, Comput. Phys. Comm. 138 (2001) 192-212]. SMMP is an open-source FORTRAN package for molecular simulation of proteins within the standard geometry model. It is designed as a simple and inexpensive tool for researchers and students to become familiar with protein simulation techniques. This announcement describes the first major revision of this software package and its newly added features.

Program summary

Title of program:SMMPCatalogue identifier:ADOJ₋v2₋0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADOJ_v2_0Program obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandOperating system under which the program has been tested:LINUX systemProgramming language used:FORTRANComputer:PC PentiumNumber of lines in distributed program, including test data, etc.:18 492Number of bytes in distributed program, including test data, etc.:278 995Distribution format:ASCIICard punching code:ASCIICatalogue Identifier of previous version:ADOJJournal Reference of previous version:F. Eisenmenger, U.H.E. Hansmann, Sh. Hayryan, C.-K. Hu, Comput. Phys. Comm. 138 (2001) 192-212Does the new version supersede the previous version?:YesNature of physical problem:Molecular mechanics computations and Monte Carlo simulation of proteinsReasons for the new version:Increased functionalitySummary of revisions:Changes in energy function and protein representation; differences in program structure and organization; new functionalities added; miscellaneous changes and additionsMethod of solution:Utilizes ECEPP2/3 and FLEX potentials. Includes Monte Carlo simulation algorithms for canonical, as well as for generalized ensemblesRestrictions on the complexity of the problem:The consumed CPU time increases with the size of protein moleculeTypical running time:Depends on the size of the molecule under simulationUnusual features of the program:No     

A fine grained parallel smooth particle mesh Ewald algorithm for biophysical simulation studies: Application to the 6-D torus QCDOC supercomputer

In order to model complex heterogeneous biophysical macrostructures with non-trivial charge distributions such as globular proteins in water, it is important to evaluate the long range forces present in these systems accurately and efficiently. The Smooth Particle Mesh Ewald summation technique (SPME) is commonly used to determine the long range part of electrostatic energy in large scale molecular simulations. While the SPME technique does not give rise to a performance bottleneck on a single processor, current implementations of SPME on massively parallel, supercomputers become problematic at large processor numbers, limiting the time and length scales that can be reached. Here, a synergistic investigation involving method improvement, parallel programming and novel architectures is employed to address this difficulty. A relatively simple modification of the SPME technique is described which gives rise to both improved accuracy and efficiency on both massively parallel and scalar computing platforms. Our fine grained parallel implementation of the modified SPME method for the novel QCDOC supercomputer with its 6D-torus architecture is then given. Numerical tests of algorithm performance on up to 1024 processors of the QCDOC machine at BNL are presented for two systems of interest, a β-hairpin solvated in explicit water, a system which consists of 1142 water molecules and a 20 residue protein for a total of 3579 atoms, and the HIV-1 protease solvated in explicit water, a system which consists of 9331 water molecules and a 198 residue protein for a total of 29508 atoms.     

NMscatt: a program for calculating inelastic scattering from large biomolecular systems using classical force-field simulations

Computational tools for normal mode analysis, which are widely used in physics and materials science problems, are designed here in a single package called NMscatt (Normal Modes & scattering) that allows arbitrarily large systems to be handled. The package allows inelastic neutron and X-ray scattering observables to be calculated, allowing comparison with experimental data produced at large scale facilities. Various simplification schemes are presented for analyzing displacement vectors, which are otherwise too complicated to understand in very large systems.

Program summary

Title of program:NMscattCatalogue identifier:ADZA_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADZA_v1_0.htmlProgram obtainable from:CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions:noNo. of lines in distributed program, including test data, etc.:573 535No. of bytes in distributed program, including test data, etc.:4 516 496Distribution format:tar.gzProgramming language:FORTRAN 77Computer:x86 PCOperating system:GNU/Linux, UNIXRAM:Depends on the system size to be simulatedWord size:32 or 64 bitsClassification:16.3External routines:LAPACKNature of problem: Normal mode analysis, phonons calculation, derivation of incoherent and coherent inelastic scattering spectra.Solution method: Full diagonalization (producing eigen-vectors and eigen-values) of dynamical matrix which is obtained from potential energy function derivation using finite difference method.Running time: About 7 hours per one k-point evaluation in sampling all modes dispersion curves for a system containing 3550 atoms in the unit cell on AMD Athlon 64 X2 Dual Core Processor 4200+.     

SMMP v. 3.0—Simulating proteins and protein interactions in Python and Fortran

We describe a revised and updated version of the program package SMMP. SMMP is an open-source FORTRAN package for molecular simulation of proteins within the standard geometry model. It is designed as a simple and inexpensive tool for researchers and students to become familiar with protein simulation techniques. SMMP 3.0 sports a revised API increasing its flexibility, an implementation of the Lund force field, multi-molecule simulations, a parallel implementation of the energy function, Python bindings, and more.

Program summary

Title of program:SMMPCatalogue identifier:ADOJ_v3_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADOJ_v3_0.htmlProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions:Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlProgramming language used:FORTRAN, PythonNo. of lines in distributed program, including test data, etc.:52 105No. of bytes in distributed program, including test data, etc.:599 150Distribution format:tar.gzComputer:Platform independentOperating system:OS independentRAM:2 MbytesClassification:3Does the new version supersede the previous version?:YesNature of problem:Molecular mechanics computations and Monte Carlo simulation of proteins.Solution method:Utilizes ECEPP2/3, FLEX, and Lund potentials. Includes Monte Carlo simulation algorithms for canonical, as well as for generalized ensembles.Reasons for new version:API changes and increased functionality.Summary of revisions:Added Lund potential; parameters used in subroutines are now passed as arguments; multi-molecule simulations; parallelized energy calculation for ECEPP; Python bindings.Restrictions:The consumed CPU time increases with the size of protein molecule.Running time:Depends on the size of the simulated molecule.     

Computational cancer cells identification by fractal dimension analysis

In the present work a software to identify cell anomaly through their fractal dimension calculation is introduced. The cell electronic microscopic image is imported to the software in gray scale and transformed to a black and white pattern in order to eliminate possible noise due to organic material close to the cell during the acquisition of the image. The number of pixels on the image contour is determined and the box-counting method is used to obtain the fractal dimension of the cell. Results for the fractal dimension have shown in very good agreement with other calculations with the advantage that the software is user-friendly, avoiding human analysis mistakes.     

REACH: A program for coarse-grained biomolecular simulation

REACH (Realistic Extension Algorithm viaCovariance Hessian) is a program package for residue-scale coarse-grained biomolecular simulation. The program calculates the force constants of a residue-scale elastic network model in single-domain proteins using the variance-covariance matrix obtained from atomistic molecular dynamics simulation. Secondary-structure dependence of the force constants is integrated. The method involves self-consistent, direct mapping of atomistic simulation results onto a coarse-grained force field in an efficient automated procedure without requiring iterative fits and avoiding system dependence.

Program summary

Program title: REACHCatalogue identifier: AEDA_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEDA_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 42 244No. of bytes in distributed program, including test data, etc.: 3 682 118Distribution format: tar.gzProgramming language: FORTRAN 77Computer: x86 PCOperating system: GNU/Linux, SUSE and Red HatRAM: Depends on the system size to be calculatedWord size: 32 or 64 bitsClassification: 3External routines: LAPACK, BLASNature of problem: A direct calculation of force field for residue-scale coarse-grained biomolecular simulation derived from atomistic molecular dynamics trajectory.Solution method: A variance-covariance matrix and the associated Hessian (second-derivative) matrix are calculated from an atomistic molecular dynamics trajectory of single-domain protein internal motion and the off-diagonal Hessian matrix is fitted to that of a residue-scale elastic network model. The resulting force constants for the residue pair interactions are expressed as model functions as a function of pairwise distance.Running time: Depends on the system size and the number of MD trajectory frames used. The test run provided with the distribution takes only a few seconds to execute.     

Simulation of distribution of random walks on a lattice

The flat histogram version of pruned and enriched Rosenbluth method (FLATPERM) is very efficient for calculating densities of classical states of polymers on a lattice. In this paper the accuracy of this method is tested by comparing simulations to an exact solution of distribution of random walks on a one-dimensional lattice. The boundary effects of restricting the region of parameter space are investigated as well. Finally by FLATPERM we calculate the distribution of end-to-end distance of self-avoiding walks on a cubic lattice and simulate the scaling behavior of most probable end-to-end distance.     

Fast clear-sky solar irradiation computation for very large digital elevation models

This paper presents a fast algorithm to compute the global clear-sky irradiation, appropriate for extended high-resolution Digital Elevation Models (DEMs). The latest equations published in the European Solar Radiation Atlas (ESRA) have been used as a starting point for the proposed model and solved using a numerical method. A new calculation reordering has been performed to (1) substantially diminish the computational requirements, and (2) to reduce dependence on both, the DEM size and the simulated period, i.e., the period during which the irradiation is calculated. All relevant parameters related to shadowing, atmospheric, and climatological factors have been considered. The computational results demonstrate that the obtained implementation is faster by many orders of magnitude than all existing advanced irradiation models while maintaining accuracy. Although this paper focuses on the clear-sky irradiation, the developed software also computes the global irradiation applying a filter that considers the clear-sky index.     

Solving Differential Matrix Riccati Equations by a piecewise-linearized method based on the conmutant equation

Differential Matrix Riccati Equations play a fundamental role in control theory, for example, in optimal control, filtering and estimation, decoupling and order reduction, etc. In this paper a piecewise-linearized method based on the conmutant equation to solve Differential Matrix Riccati Equations (DMREs) is described. This method is applied to develop two algorithms which solve these equations: one for time-varying DMREs and another for time-invariant DMREs, also MATLAB implementations of the above algorithms are developed. Since MATLAB does not have functions which solve DMREs, two algorithms based on a BDF method are also developed. All implemented algorithms have been compared, under equal conditions, at both precision and computational costs. Experimental results show the advantages of solving non-stiff DMREs and in particular stiff DMREs by the proposed algorithms.     

Efficient DNA sticker algorithms for NP-complete graph problems

Adleman's successful solution of a seven-vertex instance of the NP-complete Hamiltonian directed path problem by a DNA algorithm initiated the field of biomolecular computing. We provide DNA algorithms based on the sticker model to compute all k-cliques, independent k-sets, Hamiltonian paths, and Steiner trees with respect to a given edge or vertex set. The algorithms determine not merely the existence of a solution but yield all solutions (if any). For an undirected graph with n vertices and m edges, the running time of the algorithms is linear in n+m. For this, the sticker algorithms make use of small combinatorial input libraries instead of commonly used large libraries. The described algorithms are entirely theoretical in nature. They may become very useful in practice, when further advances in biotechnology lead to an efficient implementation of the sticker model.     

Accounting of heterogeneity in the estimation of radiation risks

The problem is considered of taking account of heterogeneity in the individual radiosensitivity in evaluating the radiation risk of the origin of a disease in a group of persons subjected to the radiation action. Mathematical models are described for the estimation of the radiation risk of the origin of the disease with and without regard for heterogeneity. The use of the described method is demonstrated with an example of the estimation of the risk of the origin of radiation-induced diseases from the class of “solid cancers” among the persons taking part in the liquidation of consequences of the accident at the Chernobyl’ Atomic Power Station (CAPS) in the years 1986–1991. It is shown that the heterogeneity disregard leads to understated estimates of the radiation risk of the origin of radiation-induced solid cancers. The effect is investigated of the length of the latent period in the origin of radiation-induced “solid cancers” in taking account of heterogeneity.     

The study of helix-coil transition of polyalanine with single histogram method

We have studied helix-coil transition in all-atom model of polyalanine in vacuum and distance-dependent dielectric function solvent model. Molecules of up to length 30 residues were investigated by Monte Carlo method. The single histogram method was performed to obtain thermodynamic quantities such as average energy and specific heat. We have also calculated probability distribution of energy and Binder cumulant as well as autocorrelation times for both case in order to identify the order of phase transition by using finite-size scaling arguments. Even though the probability distributions showed a single Gaussian peak, especially for this solvent model, the significant differences between the vacuum and solvent model were observed at their specific heat maximum values and Binder cumulant behaviors. While Binder cumulant has a deep valley for vacuum, it is not sensitive for the distance-dependent dielectric function solvent model. Our data are consistent with recently obtained results for polyalanine by multicanonical method.     

PLUMED: A portable plugin for free-energy calculations with molecular dynamics

Here we present a program aimed at free-energy calculations in molecular systems. It consists of a series of routines that can be interfaced with the most popular classical molecular dynamics (MD) codes through a simple patching procedure. This leaves the possibility for the user to exploit many different MD engines depending on the system simulated and on the computational resources available. Free-energy calculations can be performed as a function of many collective variables, with a particular focus on biological problems, and using state-of-the-art methods such as metadynamics, umbrella sampling and Jarzynski-equation based steered MD. The present software, written in ANSI-C language, can be easily interfaced with both Fortran and C/C++ codes.

Program summary

Program title: PLUMEDCatalogue identifier: AEEE_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEEE_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Lesser GPLNo. of lines in distributed program, including test data, etc.: 107 505No. of bytes in distributed program, including test data, etc.: 2 052 759Distribution format: tar.gzProgramming language: ANSI-CComputer: Any computer capable of running an executable produced by GCC compilerOperating system: Linux/UnixRAM: Depending on the number of atoms, the method chosen and the collective variables usedClassification: 23External routines: Must be interfaced with a MD code (such as GROMACS, NAMD, DL_POLY or SANDER).Nature of problem: Calculation of free-energy surfaces for biological and condensed matter systems.Solution method: Implementation of various enhanced sampling techniques.Unusual features: PLUMED is not a stand-alone program but it must be interfaced with a MD code (such as GROMACS, NAMD, DL_POLY or SANDER) that needs to be recompiled. Each interface is provided in a patch form.Running time: Depending on the number of atoms, the method chosen and the collective variables used.     

Solving Initial Value Problems for Ordinary Differential Equations by two approaches: BDF and piecewise-linearized methods

Many scientific and engineering problems are described using Ordinary Differential Equations (ODEs), where the analytic solution is unknown. Much research has been done by the scientific community on developing numerical methods which can provide an approximate solution of the original ODE. In this work, two approaches have been considered based on BDF and Piecewise-linearized Methods. The approach based on BDF methods uses a Chord-Shamanskii iteration for computing the nonlinear system which is obtained when the BDF schema is used. Two approaches based on piecewise-linearized methods have also been considered. These approaches are based on a theorem proved in this paper which allows to compute the approximate solution at each time step by means of a block-oriented method based on diagonal Padé approximations. The difference between these implementations is in using or not using the scale and squaring technique.Five algorithms based on these approaches have been developed. MATLAB and Fortran versions of the above algorithms have been developed, comparing both precision and computational costs. BLAS and LAPACK libraries have been used in Fortran implementations. In order to compare in equality of conditions all implementations, algorithms with fixed step have been considered. Four of the five case studies analyzed come from biology and chemical kinetics stiff problems. Experimental results show the advantages of the proposed algorithms, especially when they are integrating stiff problems.     

A cell multipole based domain decomposition algorithm for molecular dynamics simulation of systems of arbitrary shape

A domain decomposition algorithm for molecular dynamics simulation of atomic and molecular systems with arbitrary shape and non-periodic boundary conditions is described. The molecular dynamics program uses cell multipole method for efficient calculation of long range electrostatic interactions and a multiple time step method to facilitate bigger time steps. The system is enclosed in a cube and the cube is divided into a hierarchy of cells. The deepest level cells are assigned to processors such that each processor has contiguous cells and static load balancing is achieved by redistributing the cells so that each processor has approximately same number of atoms. The resulting domains have irregular shape and may have more than 26 neighbors. Atoms constituting bond angles and torsion angles may straddle more than two processors. An efficient strategy is devised for initial assignment and subsequent reassignment of such multiple-atom potentials to processors. At each step, computation is overlapped with communication greatly reducing the effect of communication overhead on parallel performance. The algorithm is tested on a spherical cluster of water molecules, a hexasaccharide and an enzyme both solvated by a spherical cluster of water molecules. In each case a spherical boundary containing oxygen atoms with only repulsive interactions is used to prevent evaporation of water molecules. The algorithm shows excellent parallel efficiency even for small number of cells/atoms per processor.     

TiReX: Replica-exchange molecular dynamics using Tinker

We present a driver program for performing replica-exchange molecular dynamics simulations with the Tinker package. Parallelization is based on the Message Passing Interface, with every replica assigned to a separate process. The algorithm is not communication intensive, which makes the program suitable for running even on loosely coupled cluster systems. Particular attention is paid to the practical aspects of analyzing the program output.

Program summary

Program title: TiReXCatalogue identifier: AEEK_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEEK_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 43 385No. of bytes in distributed program, including test data, etc.: 502 262Distribution format: tar.gzProgramming language: Fortran 90/95Computer: Most UNIX machinesOperating system: LinuxHas the code been vectorized or parallelized?: parallelized with MPIClassification: 16.13External routines: TINKER version 4.2 or 5.0, built as a libraryNature of problem: Replica-exchange molecular dynamics.Solution method: Each replica is assigned to a separate process; temperatures are swapped between replicas at regular time intervals.Running time: The sample run may take up to a few minutes.     

A hybrid multi-loop genetic-algorithm/simplex/spatial-grid method for locating the optimum orientation of an adsorbed protein on a solid surface

Atomistic simulation of protein adsorption on a solid surface in aqueous environment is computationally demanding, therefore the determination of preferred protein orientations on the solid surface usually serves as an initial step in simulation studies. We have developed a hybrid multi-loop genetic-algorithm/simplex/spatial-grid method to search for low adsorption-energy orientations of a protein molecule on a solid surface. In this method, the surface and the protein molecule are treated as rigid bodies, whereas the bulk fluid is represented by spatial grids. For each grid point, an effective interaction region in the surface is defined by a cutoff distance, and the possible interaction energy between an atom at the grid point and the surface is calculated and recorded in a database. In searching for the optimum position and orientation, the protein molecule is translated and rotated as a rigid body with the configuration obtained from a previous Molecular Dynamic simulation. The orientation-dependent protein-surface interaction energy is obtained using the generated database of grid energies. The hybrid search procedure consists of two interlinked loops. In the first loop A, a genetic algorithm (GA) is applied to identify promising regions for the global energy minimum and a local optimizer with the derivative-free Nelder-Mead simplex method is used to search for the lowest-energy orientation within the identified regions. In the second loop B, a new population for GA is generated and competitive solution from loop A is improved. Switching between the two loops is adaptively controlled by the use of similarity analysis. We test the method for lysozyme adsorption on a hydrophobic hydrogen-terminated silicon (110) surface in implicit water (i.e., a continuum distance-dependent dielectric constant). The results show that the hybrid search method has faster convergence and better solution accuracy compared with the conventional genetic algorithm.