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

On 7 October 2019, Prof Dr Herman Johan Christiaan Berendsen passed away shortly after his 85 birthday. Prof Berendsen is considered the founder of the molecular dynamics simulation of the biological system: the area of theoretical research that also shaped my scientific career. He worked at the University of Groningen in the picturesque Northern part of the Netherlands. There, I met him for 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 completing my doctorate, Herman allowed me to continue working in his group with a postdoc position within the “Protein Folding” EU Training network. This happened just two years before his retirement; therefore, I was also one of his last postdocs..

In 1999, when Herman retired, his colleagues at RUG and his friends organized an international conference in his honour at the University of Groningen. Recently, I could find in my historical collection of academic-related memorabilia the brochure reported below.

It was a fantastic event, and among the eminent speakers were the keynotes of two Nobel laureates: Prof. Richard Ernst and Prof. Rober Huber. After retirement, Herman dedicated himself to writing two books that distillate all his experience in molecular simulation [1] and physics 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 remembered as a great scientist and a kind and generous person. I do not remember seeing him getting irritated by anyone during a group meeting or in the famous joint BIOMOS meeting in Burg Arras (Germany) with the group of his historical collaborator Prof. Wilfred van Gunsteren, from the ETH in Zurich. He could give you sharp criticism without letting you feel like an idiot but instead try to help to solve your problem with insightful suggestions. When I was working with him in Groningen, I remember that meeting with him 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 to get good wine out of it.

Prof. van Gunsteren in an essay dedicated to his 85 birthday, recall the contribution given by Herman to the establishment of the computer simulation of the biological system [3]. It 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 a biological system.

Since then, Prof. Berendsen has provided fundamental contributions to developing the molecular dynamics simulation of biological molecules method. For this reason, he is in my personal opinion one of the founder of the MD simulation of biological systems. Here, I will shortly recall and describe some of Herman’s most important contributions to the field.

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 pressure of the system to reproduce the same experimental conditions. In the 80s of the last century, Berendsen invented and simple and robust algorithm to achieve this purpose. He called this method the weak-coupling method to an external thermal bath, but it soon became known as the Berendsen method. So far, it is one of the most widely used 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 a 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 the 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 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) to get 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 Physics23 (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 integration step length is limited by the frequency of the faster oscillator present in the molecule. This frequency is related to the type of vibration and the oscillator’s mass. In fact, to integrate these motions correctly (assuring the conservation of the total energy of the system), it is necessary to choose values of the integration step that are very small. If we are not interested in studying the vibrational properties of the system, then it is an excellent choice to remove these vibrations using holonomic constraints. Using these constraints has the advantage of allowing larger time steps (up to 10 times), speeding up the simulation. At the dawn of the MD simulation of biomolecules, as it is nowadays, extending the simulation time allows the study of dynamic processes occurring at different time scales. Berendsen and collaborators developed the SHAKE algorithm to constrain 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 are varied until all the constraint conditions are satisfied, and a position \vec{r}_i(t+\Delta t) is found. This procedure is performed by correcting the constrained atoms’ velocities after applying the SHAKE algorithm to eliminate velocity components along the direction of the constraint. Using 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. 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 the 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 Chem26(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 the beginning of 1990, Prof Berendsen initiated the development at the University of Groningen 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 [4]. The specific parallel program for MD simulations written in the C programming language was also written by David van der Spoel, a talented PhD student now a full professor at Uppsala University in Sweden. Unfortunately, the rapid evolution of low-cost parallel computing based on PC clusters made it disadvantageous to continue developing the specialized Gromacs hardware. Nevertheless, the Gromacs software package was initially designed and implemented on different parallel architectures based on Silicon Graphics clusters. The program’s first release was the 1.6, released on 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 development of the Linux OS on ordinary PCs rapidly spread the program’s use in the scientific molecular simulation community. The introduction of new intelligent ideas for calculating the non-bonded pairs list, the virial of the forces and using assembler-level optimized routines for the non-bonded interaction calculations on Intel and AMD processors enormously improved the program’s performance on affordable PC clusters and, consequently, its extensive use. In the spirit of the freeware movement, Herman wanted that the program remains public domain and free and available to everyone working in academics. Since 2001, the program has been licensed as GNU Public Licence (GPL) freeware, and its source code is still freely available and open for development. David van der Spoel and his collaborators continue to maintain and develop it at Uppsala University with other contributors worldwide. As stated in the current documentation, the name Gromacs is not associated 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 would be too long to report here. I have co-authored with him 5 papers [5-9] and various poster abstracts [10-14], and I am delighted that one of them is also my first published one [5].

I got to know Herman at the very end of his carrier. So this is just a personal and limited tribute to his intellectual achievements in the field. His PhD students and early co-authors can provide a better picture of him as a supervisor and scientist. However, there is no doubt that his legacy will continue to impact the field of MD simulations for many years.

After I wrote the first version of this blog, Macuglia, Roux, and Ciccotti published a comprehensive historical review where the interested reader can find a more extensive and first person description of all actors that have played the development of the MD simulation of biological molecules[15].


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 BioMolecular Simulation: The EightWeek CECAM Workshop ‘Models for Protein Dynamics’ of 1976. Helvetica Chimica Acta 102(8), 2019.
  4. Bekker, H., Berendsen, H.J.C., Dijkstra, E.J., Achterop, S., van Drunen, R., Van der Spoel, D., Sijbers, A., Keegstra, H., Reitsma, B. and Renardus, M.K.R., 1993. Gromacs: A parallel computer for molecular dynamics simulations. In Physics computing (Vol. 92, pp. 252-256).
  5. A. Di Nola, D. Roccatano, H. J. C. BerendsenMolecular dynamics simulation of the docking of substrates to proteins. Proteins: Struct., Funct., Genet., 19, 174-182 (1994).
  6. A. Amadei, D. Roccatano, M. E. F. Apol, H. J. C. Berendsen, A. Di Nola. Prediction of the liquid-vapor equilibrium pressure using the Quasi-Gaussian entropy theory. J. Chem. Phys., 105, 7022-7025 (1996).
  7. D. Roccatano, H. J. C. Berendsen, P. D’Angelo. Assessment of the validity of intermolecular potential models used in molecular dynamics simulations by extended X-ray absorption fine structure spectroscopy: A case study of Sr in methanol solution. J. Chem. Phys., 108, 9487-9497 (1998).
  8. D. Roccatano, A. Amadei, M. E. F. Apol, A. Di Nola, H. J. C. Berendsen. Application of the quasi-Gaussian entropy theory to molecular dynamics simulations of Lennard-Jones fluids. J. Chem. Phys., 109, 6358-6363 (1998).
  9. D. Roccatano, A. Amadei, A. Di Nola, H. J. C. Berendsen. A molecular dynamics study of the 41-56 β-hairpin from B1 domain of Protein G. Protein Sci., 8, 2130-2143 (1999).
  10. D. Roccatano, A. Amadei, A. Di Nola, H. J. C Berendsen. Study of a β-hairpin forming peptide in aqueous solution using molecular dynamics simulations. Journal of Bioscience, 24/S1, 49 (1999). POSTER abstract.
  11. D. Roccatano, A. Amadei, H.J.C Berendsen. Study of Peptide Unfolding Processes using Molecular Dynamics simulation and Essential Dynamics Analysis.EMBO Workshop. Protein Folding & Misfolding inside and outside the Cell. St, Catherine’s College, University of Oxford, UK, March 25-28th, 1998. POSTER abstract.
  12. D. Roccatano, H.J.C Berendsen. A Molecular Dynamics Simulation Study of Interleukin 1β receptor. Lunteren Meeting. Dec 7-8, 1998. Lunteren (The Netherland). POSTER abstract.
  13. D. Roccatano, A. Amadei, H.J.C Berendsen. Study of Peptide Unfolding using Molecular Dynamics simulation and Essential Dynamics Analysis. SON-working party for Protein Research. Lunteren Meeting. Dec 8-9, 1997. Lunteren (The Netherland). POSTER abstract.
  14. D. Roccatano, A. Amadei, H. J. C Berendsen. Theoretical approach to the study of local density fluctuations in simulated fluids. Stichting Scheikunding Onderzoek in Nederland. Wetenschappelijke bijeenkomst Werkgemeenschap voor Vloeistoffen en Grensvlakken. 17-18 Febraury 1997, Lunteren (The Netherland). POSTER abstract.
  15. Macuglia, D., Roux, B. and Ciccotti, G., 2022. The emergence of protein dynamics simulations: how computational statistical mechanics met biochemistry. The European Physical Journal H47(1), 13.

2 thoughts on “A personal tribute to one of the founders of MD simulation of biological molecules: Prof Herman J.C. Berendsen (1934-2019)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.