A personal tribute to the founder of MD simulation of biological molecules: Prof Herman J.C. Berendsen (1934-2019)

Posted on October 23, 2019 by

On the 7 October 2019, Prof Dr Herman Johan Christiaan Berendsen passed away just shortly after his 85 birthday. Prof Berendsen is considered the founder of the molecular dynamics simulation of biological system: the area of theoretical research that also shaped my scientific career. He was working at the University of Groningen in the pictoresque northen part of the Netherlands. It was there that I meet him the first time as it allowed me to conduct research in his lab during the last year of my doctorate researches training at the University of Rome “La Sapienza”. After I completed my doctorate, Herman gave me the opportunity to continue working in his group with a postdoc position within the “Protein Folding” EU Training network. This happen just two years before his retirement and therefore I was also one of his last postdocs. After retirement, Herman dedicated himself to write two books that distillate all his experience in the area of molecular simulation [1] and in the education [2]. He stated in a project on the social scientific network Researchgate that “I am retired and work occasionally on methods for multiscale simulations.”

Prof. Berendsen will be remember as a great scientist and a kind and generous person. I do not remember seen him getting irritated by anyone during group meeting, or in the famous joint BIOMOS meeting in Burg Arras (Germany) with the group of his historical collaborator Prof. Wilfred van Gunstern from the ETH in Zurich. He could give you sharp criticism without let you feel an idiot but instead trying to help to solve your problem with insightfully suggestions. When I was working with him in Groningen, I remember that having meeting with him it was like going to have a good harvest of ideas: he was always ready to give you the best grapes but it was up to you get good wine out of it.

Prof. van Gunsteren in a essay dedicated to his 85 birthday recall the contribution given by Herman to the establishment of the computer simulation of biological system [3]. All started with the CECAM conference in Paris organized by Herman in 1976, gathering together 20 scientists, two of them awarded in 2013 with the Nobel prize in Chemistry for their contribution to the field of computer simulations of biological system.

Since then Prof. Berendsen has given many contribution to the development of the molecular dynamics simulation of biological molecules. Here, I am just going to recall and describe some of the most important.

Berendsen Thermostat and Barostat

H. J. C. Berendsen, J. P. M. Postma, A. DiNola and J. R. Haak. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81 (1984) pp. 3684-3690.

According to Web of Science (accessed on Nov 2019), this is the most cited paper of Herman with 18242 citations.

In molecular dynamics simulation, you need to control the temperature and the pressure of the system to reproduce the same experimental conditions. In the 80s of the last century, Berendsen invented and simple and robust algotothm to achieve this purpose. He called this method the weak-coupling method to an external thermal bath but it soon become known ad the Berendsen method. So far, it is one of the most widely used  in  MD simulations. One of the advantages of this method is its simplicity. The algorithm simulates the coupling of the system to an external thermal bath with fixed temperature. The equations of the motion are modified in such a way that the current temperature of  the system T has a first order relaxation  towards the reference value T_0: 

\frac{dT}{dt} =  \frac{1}{\tau} [T_0-T(t)].

The intensity of the interaction between the thermostat bath and the system is determined by the time constant, \tau. Schematically the coupling is obtained by modifying the velocity of the atoms at each step.  This can be done by multiplying the velocities of  each particle by a constant \lambda. The new temperature is given by the relation

T'(t)= \sqrt{\left[1+  \frac{\Delta t}{\tau}\left(\frac{T_0}{T(t)}-1\right) \right]} T(t)

The weak-coupling method is also used in order to maintain to pressure of the system at a constant value. The equations of the motion are modified in such a way that the pressure of  the system has a first order relaxation  towards the reference value: 

\frac{dP(t)}{dt} = \frac{P_0-P(t)}{\tau_p}

the hydrostatic pressure P is defined for an isotropic system by the relation:

P(t) = \frac{2}{3}\frac{E_{K}(t) - \Theta(t)}{V_{s}(t)}

in which V_{s}  represents the volume of the system and the virial of intermoleculat forces \theta(t). The isothermal compressibility (\beta_T) is  the proportionality constant that relates the change of pressure  to the change of volume by the relation:

\Delta P(t) = -\frac{\Delta V(t)}{\beta_T V(t)}.

The change of volume can be obtained by multiplying it by a factor \mu:

\Delta V(t) = (\mu(t)^3 - 1 )V(t).

given by the expression

\mu(t) = \left\{ 1- \beta_T \frac{\delta t}{\tau_P} \left[ P(t)-P_0 \right] \right\}^{1/3}.

Practically this is accomplished by rescaling the coordinate of the atoms in the system, \vec{r}_i(t), using the same value of \mu

\vec{r}_i^N(t)=\mu \vec{r}_i(t)

According to the equation, the relaxation of the pressure $P(t)$ towards the reference value is controlled by the relationship between the isothermal compressibility of the system and the value chosen for relaxation time of the pressure (\tau_p).  It is not necessary to know the value of \beta_t exactly, since \tau_p is a parameter that can always be adjusted to obtain the proper reference pressure. Such a value should be sufficiently small (strong coupling) in order to obtain the equilibrium value of the pressure, but large enough (weak coupling) to avoid perturbing the properties of the system.

Development of the SHAKE algorithm

Ryckaert, J-P; Ciccotti G; Berendsen HJC (1977). “Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes”. Journal of Computational Physics. 23 (3): 327–341.

According to Web of Science (accessed on Nov 2019), this is the second most cited paper of Herman with 13482 citations.

For flexible molecules the length of the  integration step is limited by the frequency of the faster oscillator present in the molecule. This frequency is related to the type of vibration and to the mass of the oscillator. In fact, in order to integrate these motions correctly (assuring the conservation of the total  energy of the system), it is necessary to choice values of the  integration step that are very small. In the case, we are not interested to study vibrational properties of the system, than it is a good choice to remove these vibration using  holonomic constraints. The use of these constraints has the advantage to allow the use of larger time steps (till to 10 times), speading up the simulation. Nowadays as (and even more) at the dawn of MD simulation of biomolecules, extend the simulation time span allow to access the study of dynamic processes a occurring at different time scale. Berendsen and collaborator developed an algorithm, named SHAKE, to constraints bond distances based on the coordinate relaxation method. The SHAKE is an iterative procedure in which after each integration step, the positions \vec{r}_i(t+\Delta t) of atoms is varied until all the constraint conditions are satisfied, and a positions \vec{r}_i(t+\Delta t) is found. This procedure is  performed correcting the velocities of the  constrained  atoms, after the application of the SHAKE algorithm, to eliminate velocity  components along the direction of the constraint.  The use of the SHAKE algorithm allows increasing the integration step in MD simulations of $2-4$ times.

Simple Point Charge (SPC) models of water

H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, The missing term in effective pair potentials, Journal of Physical Chemistry 91 (1987) 6269-6271.

According to Web of Science (accessed on Nov 2019), this is the most cited paper of Herman with 8023 citations.

Realistic simulations of biomolecules require the presence of an explicitly solvent model. The model of water parameterized by Berendsen and his group in Groningen is a simple rigid model. Still, capable of reproducing with reasonable accuracy the experimental density, enthalpy of evaporation and the variant SPC/E also the diffusion coefficient of water at 298 K. Given its simplicity, the models are still broadly used for simulation of biological molecules. The parameters of the model are reported in the following figure.

Figure 1: SPC and SPC/E water models

The GROMACS Package for MD Simulations

H. J. C. Berendsen, D. van der Spoel and R. van Drunen. GROMACS: A message-passing parallel molecular dynamics implementation. Comp. Phys. Comm. 91 (1995) pp. 43-56.

Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark A.E, Berendsen H.J.C (2005). “GROMACS: fast, flexible, and free”. J Comput Chem. 26(16): 1701–18

According to Web of Science (accessed on Nov 2019), these two papers are the sixth and fifth most cited with 5048 and 7048 citations, respectively. 

The name of this program is an acronym for Groningen Machine for Chemical Simulations. At beginning of 1990, Prof Berendsen initiated the development at University of Groningen (Netherlands) of a parallel computer called Groningen Machine for Chemical Simulations or Gromacs ring processor. It was a parallel computer based on the i860 transputers and specifically designed to perform MD simulations with high computational efficiency. The specific parallel program for MD simulations written in the C programming language was also written by mainly David van der Spoel a talented PhD student now full professor at the University of Stockhom in Sweden. The rapid evolution of low-cost parallel computing based on PC clusters made disadvantageous to keep continuing developing the specialized Gromacs hardware. Neverthenless, the Gromacs software package was further developed and implemented on different parallel architectures initally based on Silicon Graphics clusters. The first release of the program was the 1.6 and it was released on the 27/11/1997 (see in Figure 1 the original webpage and in Figure 2 an example of header comments in the source codes).

Figure 2: the first web page of the Gromacs 1.6 distribution.
Figure 3: Header comments in the source code file mdrun.c of Gromacs 1.6.

The contemporary developed of the Linux OS on ordinary PC rapidly spread the use of the program in the scientific community of molecular simulation. The introduction of new smart ideas for the calculation of the non-bonded pairs list, the virial of the forces and the use of assembler level optimized routines for the non-bonded interaction calculations on Intel and AMD processors enormously improved the performance of the program on affordable PC clusters and consequently its extensive use. In the spirit of the freeware movement of Herman wanted that the program remains public domain free available to everyone working in academics. Since 2001, the program has been licensed as GNU General Public Licence (GPL) freeware and its source code is still freely available and open for development. Currently, David van der Spoel and his collaborators continue to maintain and develop it at the University of Upsala with other contributors worldwide. As stated in the current documentation, the name Gromacs has not associated anymore with its original acronym, another example of Herman’s scientific generosity and modesty. The official website of Gromacs is http://www.gromacs.org. 

Herman, with his group, has contributed to many other breakthrough achievements in MD simulation of a biological system that it would be too long to report here. I have co-authored 5 papers with Herman, and I am delighted that one of them is also my first published one [4].

There is no doubt that his legacy will continue to impact the field of MD simulations for many years to come.

REFERENCES

  1. H.J.C. Berendsen. Simulating the Physical World: Hierarchical Modeling from Quantum Mechanics to Fluid Dynamics. Cambridge University Press, (12 July 2007).
  2. H.J.C. Berendsen.A Student’s Guide to Data and Error Analysis (Student’s Guides).Cambridge University Press, (7 April 2011).
  3. W. van Gunsteren. The Roots of Bio‐Molecular Simulation: The Eight‐Week CECAM Workshop ‘Models for Protein Dynamics’ of 1976. Helvetica Chimica Acta 102(8), 2019.
  4. A. Di Nola, D. Roccatano, H. J. C. Berendsen. Molecular dynamics simulation of the docking of substrates to proteins. Proteins: Struct., Funct., Genet., 19, 174-182 (1994).

Posted in Research, What is new | Leave a comment

The First 150 Years of the Periodic Table of the Elements

Posted on September 1, 2019 by

That the nobility of man, acquired in a hundred centuries of trial and error, lay in making himself the conquerer of matter, and that I had enrolled in chemistry because I wanted to remain faithful to this nobility. That conquering matter is to understand it, and understanding matter is necessary to understanding the universe and ourselves: and that therefore Mendeleev’s Periodic Table, which just during those weeks we were laboriously learning to unravel, was poetry, loftier and more solemn than all the poetry we had swallowed down in liceo; and come to think of it, it even rhymed!

Primo Levi, The Periodic Table.

This year marks the 150th anniversary of the periodic table of the elements (TPE) which currently has 118 entries, the latest arrival (the Tennessium) was discovered 10 years ago (2009), and I feel obliged as a chemist to give some a small informative contribution to celebrate this important event.

Continue reading →
Posted in Uncategorized | Leave a comment

Retro programming nostalgia III: the MSX Microcomputer and the Orbit of the Planets in the Solar System

Posted on August 31, 2019 by

In a recent article, I have explained the Euler’s method for solving ordinary differential equationsusing as a motivation the fictionalized version in the film Hidden Figures of the scientific contribution of Katherine Goble and her two colleagues to the NASA space program. I have also described as an example of application a program in awk programming language for calculating the orbits of planets of the solar system. My interest in astrodynamics come back to my juvenile age, when still going to high school, my parents decided to gift me a more sophisticated microcomputer than my previous one (the celebrated Commodore VIC 20). So I became a programmer of a Philips MSX VG 8010 that I still jelously own in its original box. So, powered by the versatile Federico Faggin’s Zilog Z80 processor with a clock 3.58 MHz, with an impressive (for a previous owner of a VIC20 with a mere 3.583 kB!) memory of 32 kB RAM , 16kB of video RAM and a dedicated tape-record device as storage system, I started to write more sophisticated in MSX Basic. At that time, I was eagerly following the department “Ricreazioni al Computer” by the famous computer scientist A. K. Dewdney on the magazine “Le Scienze”, the Italian edition of Scientific American. The new microcomputer allowed me to experiment with the fascinating computational topics that Dewdney was offering every month. One of these topics was dedicated to the simulation of stars using the algorithm based on the Euler integration of the Newton equation. Following the instruction of Dewney, I managed to write a small program in MSX basic and this was the starting of my interest in computational astronomy.

Continue reading →
Posted in Programming, What is new | Leave a comment

Numerical Integration of Differential Equations. Part I.: Katherine Goble and the Euler’s Method.

Posted on May 31, 2019 by

This blog was inspired by the beautiful 2016 movie Hidden Figures (based on the book of the same name by M. L. Shetterley) which tell the dramatic story of three talented black women scientist that worked as “human computers” for NASA in 1961 for the Mercury project.

Figure 1: Official theatrical poster of the movie and the phFoto of the real protagonist. From left to right. Mary Jackson, Katherine Goble and Dorothy Vaughan. (source: wikipedia)

In the movie, the mathematician Katherine Goble (interpreted by Taraji P. Henson), had a brilliant intuition on how to numerically solve the complex problem to find the transfer trajectory for the reentry into the Earth atmosphere of the Friendship 7 capsule with the astronaut John Glenn on board. In the particular scene, she was standing together with other engineers and the director of the Langley Research Center (a fictional character interpreted by Kevin Coster) in front of the vast blackboard looking to graph and equations when she says that the solution might be in the “old math” and she runs to take an old book from a bookshelf with the description of the Euler method. The scene is also nicely described in the youtube video lesson by Prof. Alan Garfinkel of the UCLA. A detailed description of the numerical solution based on the original derivation of K. Globe is in the Wolfram blog website.

Katherine Globe was using for these complex calculation her brilliant brain with the support of a mechanical calculator (the Friden STW-10, in the movie, this machine is visible in different scenes). In a scene of the film, she revealed that her typical computing performance was of 10000 calculations per day and probably for calculations, she was not referring to single arithmetic operations! These exceptional mathematical skills have given a significative contribution at the beginning of the American space program, but it became insufficient to handle the more complex mathematics necessary to land the man on the Moon, and the other fantastic NASA achievements.

Continue reading →
Posted in Programming, Science Topics, What is new | Leave a comment

The Calculation of the Lattice Energy: The Born-Haber Cycle

Posted on May 26, 2019 by

My blog in italian on this topics is very popular and for this reason I decided to add an English translation (when I have some free time, I will also translate the text in the Figure and Table). So be tune and more will come!

The stability of a crystal lattice at constant T and P conditions is linked to the Gibbs free energy of lattice formation by the relations

M^+ (g) + X^- (g) \rightarrow MX (s) \hfill (1)

\displaystyle \Delta G^{\circ} = \Delta H^{\circ} - T \Delta S^{\circ} \hfill (2)

If {\Delta G^{\circ}} is more negative for the formation of the {I} structure than for the {II} structure, the {II \rightarrow I} transition will be spontaneous and the solid will have that structure.

Continue reading →
Posted in Science Topics, What is new | Leave a comment

Modelling Natural Patterns and Forms II: (Easter) Eggs

Posted on April 2, 2019 by

When you start with a portrait and search for a pure form, 
a clear volume, through successive eliminations, you arrive 
inevitably at the egg. Likewise, starting with the egg and 
following the same process in reverse, one finishes with 
the portrait.

PABLO PICASSO

Easter is coming and what better time to talk about eggs!

During my recent mathematical explorations of natural shapes and forms, my attention has been catched by the shape of birds eggs. In the interesting book by J. Adams, A Mathematical walk in Nature [1], you can find a short review on the different mathematical modelling approach to describe the shape of an egg. Among them, the geometrical one by Baker [2] is revealed one of the most versatile as it can very accurately reproduce the shapes of a large variety of bird eggs [2]. More recently, the model was used to perform a systematic and comparative study of the shape of bird eggs. This study, published on Science magazine [3], a two-dimensional morphological space defined by the parameters of the Baker’s equation, has been used to show the diversity of the shape of 1400 species of birds. Combining these information with a mechanical model and phylogenetics information, the authors have shown that egg shape correlates with flight ability on broad taxonomic scales. They concluded that adaptations for flight may have been critical drivers of egg-shape variation in birds [3].

Continue reading →
Posted in Leonardo's Corner, Science Topics, What is new | Leave a comment

I Primi 150 Anni della Tavola Periodica degli Elementi

Posted on January 28, 2019 by

Che la nobiltà dell’Uomo, acquisita in cento secoli di prove e di errori, era consistita nel farsi signore della materia, e che io mi ero iscritto a Chimica perché a questa nobiltà mi volevo mantenere fedele. Che vincere la materia è comprenderla, e comprendere la materia è necessario per comprendere l’universo e noi stessi: e che quindi il Sistema Periodico di Mendeleev, che proprio in quelle settimane imparavamo laboriosamente a dipanare, era una poesia, più alta e più solenne di tutte le poesie digerite in liceo: a pensarci bene, aveva perfino le rime! 

Primo Levi, Il sistema periodico

Il 6 Marzo del 1869 il chimico russo Dmitri Ivanovich Mendeleyev presento’ alla Societa’ di Chimica Russa, una comunicazione dal titolo La dipendenza delle proprieta’ degli elementi chimica dal peso atomico. In questa storica comunicazione, Mendeleev pesento’ una tabella in cui organizzava gli elementi chimici allora noti. Questa tabella segno’ anche la fama del suo autore poiche’ fu la prima versione della moderna tavola periodica degli elementi chimici.

La tavola periodic degli elementi di Mendeleev. I trattini rappresentano elelmenti sconosciuti nel 1871. (fonte della figura: wikipedia)

Mendeleyev, preparando una seconda edizione del suo libro di chimica, stava cercando un modo per classificare gli elementi chimici allora conosciuti (53 ovvero meno della meta’ di quelli che conosciamo oggi) per fare chiarezza sulle loro proprieta’. In una nota, Mendeleyev racconta che l’ispirazione gli sia venuta in sogno (non e’ la prima volta che Orfeo suggerisce a chimici le loro grandi scoperte scientifici!) [2]:

I saw in a dream a table where all the elements fell into place as required. Awakening, I immediately wrote it down on a piece of paper.

Continue reading →
Posted in Leonardo's Corner, Science Topics, What is new | Tagged , , | Leave a comment