The emergence of protein dynamics simulations: how computational statistical mechanics met biochemistry

Abstract

In this essay, we aim to illustrate how Martin Karplus and his research group effectively set in motion the engine of molecular dynamics (MD) simulations of biomolecules. This process saw its prodromes between 1969 and the early 1970s with Karplus’ landing in biology, a transition that came to fruition with the treatment of 11-cis-retinal photoisomerization and the development of an allosteric model to account for the mechanism of cooperativity in hemoglobin. In 1977, J. Andrew McCammon, Bruce Gelin, and Martin Karplus published an article in Nature reporting the MD simulation of bovine pancreatic trypsin inhibitor (BPTI). This publication helped initiate the merger of computational statistical mechanics and biochemistry, a process that Karplus undertook at a later stage and whose beginnings we propose to reconstruct in this article through unpublished accounts of the key people who participated in this endeavor.

1 Introduction

The strand of biological research that Martin Karplus (born 1930) began pursuing in the early 1970s culminated in one of the first high-impact applications of contemporary computer technologies to the life sciences. Through the innovative use of cutting-edge computers for that time (e.g., IBM 7090, IBM System/370 Model 168), his group applied techniques for the deployment of ab initio and semiempirical quantum mechanics, as well as classical dynamics and statistical mechanics, to study chemical processes pertaining to biological systems, in essence, turning the focus of this novel computational approach resolutely toward biochemistry.¹ Indeed, the re-emergence of classical physics in a world of theoretical chemistry then dominated by quantum mechanics represents a rather peculiar circumstance.² In this article we aim to explore specific aspects of this occurrence, pondering how physical science, computer technology, and biology came together in Karplus’ scientific pursuits, ushering in innovative approaches to the life science.

The journey that led Karplus and his group to the molecular dynamics (MD) of biomolecules began at the Weizmann Institute of Science, Israel, between 1969 and 1970. This time was a period of exploration in which Karplus, essentially new to biology, kept abreast of the state-of-the-art in calculating the ground state of simple molecular aggregates from the Born–Oppenheimer (BO) perspective.³ This was followed by a second phase, begun as early as 1971 and essentially terminated in 1975, in which his laboratory applied itself to the study of the geometric structure and properties of 11-cis-retinal (retinaldehyde), the chromophore of the visual pigment, and to the development of a statistical mechanical model for the cooperativity of hemoglobin. Finally, there was the beginning of MD studies of biomolecules, culminating in 1977 with the publication in Nature of the bovine pancreatic trypsin inhibitor (BPTI) simulation, consecrating an area of research in which Karplus would be involved for the rest of his scientific career.⁴ Karplus and his group were able to show that proteins are flexible and move continuously and that their internal structure is characterized by statistical fluctuations at room temperature. Physicists already had experience understanding this type of behavior via molecular dynamics simulations of liquids, and their knowledge was rapidly transferred to biomolecular simulations. It was in this way that statistical physics made its grand entry into the domain of biochemistry, as biochemists realized that understanding biomolecular dynamics could only be addressed reliably and manageably with computational statistical mechanics. In addition to guiding the reader through an analysis of the above periods, we aim to highlight a pattern that was repeated in each of them. On the one hand, there was the persistent desire to extend to increasingly complex molecular systems the scheme of fundamental computer simulations (FCS) that he had co-developed in the 1960s to study the scattering reaction (H, ${\rm H}_{2}$) and, before that, the potential energy surface for ${\rm H}_{3}$, its unstable intermediate state.⁵ On the other, there was the willingness and ability to connect with young researchers of great acumen, whom Karplus enlisted in his group at Harvard University and put to work on the programs mentioned above. Our contribution includes some direct interactions with Martin Karplus and the people who participated in his venture, whose testimonies are detailed and commented on below.

2 The background

When we talk about molecular simulations (MS) we refer to an atomic approach that allows us to treat the dynamical evolution of a given atomic or molecular system to derive some of its relevant microscopic and macroscopic properties.⁶ When applied to biological macromolecules, MS provide an approach to the motions of the atoms composing complex macromolecular structures like proteins or DNA, generally in solution, yielding insights “into biological phenomena such as the role of flexibility in ligand binding and the rapid solvation of the electron transfer state in photosynthesis.”⁷ MS are part of a broader category of simulations, which are called computer simulations, and our interest lies in the MS techniques that belong to FCS, i.e., those that start from the fundamental laws of physics and proceed by solving the equations that describe the time-evolution of complex dynamical systems.⁸ In general, these equations cannot be solved using a paper-and-pencil approach, i.e., analytically, but require advanced computations.⁹ From this perspective, computers allow us to obtain numerical solutions, letting us “see” mathematical models in action, and make predictions about the behavior of the systems under investigation.¹⁰

Overall, the field of MS is subsumed into the two main offshoots represented by Metropolis Monte Carlo (Metropolis MC) and molecular dynamics (MD), both of which constitute classes of computational approaches through which scientists can achieve the goals of MS.¹¹ The MC methods are based on the idea that “a determinate mathematical problem can be treated by finding a probabilistic analogue which is then solved by a stochastic sampling experiment.”¹² Such an experiment involves “the generation of pseudo-random numbers followed by a few arithmetic and logical operations, which are often the same at each step.”¹³ This approach, devised in 1946 by Stanisław Ulam (1909–84) and developed in collaboration with John von Neumann (1903–57) in 1947 using the ENIAC, was enhanced considerably in the early 1950s, with the development by Marshall Rosenbluth (1927–2003) et al. of the Metropolis MC method.¹⁴ This happened at Los Alamos through the use of the MANIAC computer, of whose construction Nicholas Metropolis (1915–99) was the director.¹⁵ In contrast, MD involves “the solution of the classical equations of motion (Newton’s equations) for a set of molecules.”¹⁶ This task was first accomplished, for a hard-sphere system, by Berni Alder (1925–2020) and Thomas Wainwright (1927–2007) using an IBM 704 at the University of California Radiation Laboratory at Livermore. The two presented their results at an international symposium on transport processes in statistical mechanics organized in 1956 in Brussels by Ilya Prigogine (1917–2003) and published them in 1957.¹⁷ The MD method was greatly extended by Aneesur Rahman (1927–87) in 1964 who worked on particles characterized by continuous potentials, succeeding in making the major contribution that led to the simulations of realistic atomic models.¹⁸ Rahman was then joined by Frank Stillinger (born 1934) to implement the first MD simulation of liquid water.¹⁹ An important aspect of the field of MS is its deep connection with statistical mechanics: in a sense, MS represent computer-implemented “brute force” statistical mechanics. In particular, MD simulations are:

[…] the computer approach to statistical mechanics [italics added]. As a counterpart to experiment, MD simulations are used to estimate equilibrium and dynamic properties of complex systems that cannot be calculated analytically. […] MD simulations occupy a venerable position at the crossroads of mathematics, biology, chemistry, physics, and computer science.²⁰

While a narrow class of problems in statistical mechanics with idealized Hamiltonian (e.g., the Einstein crystal, the Ising model in one and two dimensions) are exactly solvable in the sense that one can achieve an accurate specification of the observable macroscopic properties in closed-form, the vast majority of problems—perhaps the most interesting ones—are actually not.²¹ Metropolis MC techniques show a robust ability to analyze thermodynamic equilibrium but not dynamical phenomena. On the other hand, MD methods can be used for tackling the properties of systems in both equilibrium and non-equilibrium, i.e., dynamic, situations. The two techniques differ in their assumptions and range of applicability, although in the case of equilibrium conditions, Metropolis MC methods can be in some cases much faster.

3 The Weizmann Institute of Science

In the fall of 1966, Karplus became professor at Harvard University and focused his attention mainly on scattering theory and many-body perturbation theory (MBPT).²² Until the early 1970s, his research had a predominance of electronic-structure quantum mechanical subjects, e.g., nuclear magnetic resonance (NMR), electron paramagnetic resonance, along with a major component on reaction kinetics.²³ At Harvard, he also worked on reaction dynamics and hyperfine interactions in electron spin resonance, but apparently it was not easy for him to move away from pure theoretical chemistry.²⁴ At some point, however, he realized that his approach needed to change and the time became ripe to make the transition: “After I had been at Harvard for only a short time, I realized that if I was ever to return to my long-standing interest in biology I had to make a break with what had been thus far a successful and very busy research program in theoretical chemistry.”²⁵

It was July 1969, when he embarked on a six-month stay at the Weizmann Institute of Science in Rehovot, Israel, a leading research institute promoting interdisciplinary research in the natural and exact sciences.²⁶ Karplus chose Israel because of the well-stocked Weizmann Institute library as well as the possibility of establishing a connection with chemical physicist Shneior Lifson (1914–2001), who at the time was working on polymer theory and the development of empirical potential energy functions to describe the energy surfaces of small molecules.²⁷ Karplus was looking for a way to connect his expertise with biology, and take a moment of pause to reflect, just as he did at the University of Oxford when he worked under Charles Coulson (1910–74) as a National Science Foundation (NSF) postdoctoral fellow from 1953 to 1955.²⁸ “The sabbatical [at the Weizmann Institute],” Karplus remarked, “gave me the possibility […] to read and explore a number of areas in which I hoped to do constructive research by applying my expertise in theoretical chemistry to biology. Discussions with Lifson and many visitors to his group helped me in these explorations.”²⁹ There was one volume, in particular, that proved to be of great stimulus to him at that time. The book, Structural Chemistry and Molecular Biology, written to celebrate Linus Pauling’s (1901–94) sixty-fifth birthday, contains an essay by Ruth Hubbard (1924–2016) and George Wald (1906–97) entitled “Pauling and the Carotenoid Stereochemistry.” There is one passage in it, in particular, that places great emphasis on the markedly quantitative approach that Pauling adopted in his research, and which ended up greatly influencing Karplus, as we will show later in our analysis:

One of the admirable things about Linus Pauling’s thinking is that he pursues it always to the level of numbers. As a result, there is usually no doubt of exactly what he means. Sometimes his initial thought is tentative because the data are not yet adequate, and then it may require some later elaboration or revision. But it is frequently he who refines the first formulation.³⁰

The Weizmann Institute was an attractive intellectual arena that eventually influenced the rest of Karplus’ career, if only because of two key people who would be directly involved in his biological research: one was Barry Honig (born 1941), the other was Arieh Warshel (born 1940).³¹ The former had received his doctorate from the Weizmann Institute in 1968, working under Tel Aviv University physical chemist Joshua Jortner (born 1933).³² In conversation with the authors, Honig recounted that despite his background in physical chemistry, he felt a strong desire to devote himself to biology—a trend followed by other scientists who had been around Karplus in the early stages of his career, such as Max Delbrück (1906–81) and Seymour Benzer (1921–2007), who were both trained in physics.³³ Honig’s words provide a direct insight into the spirit of enthusiasm surrounding biology at the time and the fascination it held for scientists educated in other disciplines, particularly physics. In Israel, Honig, who had been working on spectroscopy as the topic of his doctoral dissertation, applied for an NSF fellowship to focus on energy transfer in photosynthetic processes. He got the fellowship, and, in the summer of 1968, he moved to Harvard to work with Karplus.³⁴ It was there that the latter directed Honig’s attention to retinal. Before leaving for Israel, in fact, Karplus suggested that Honig read the article by Hubbard and Wald which appeared in Structural Chemistry and Molecular Biology.³⁵ Warshel, for his part, was a graduate student who obtained his Ph.D. in chemical physics under Lifson in 1969. He had already collaborated with biophysicist Michael Levitt (born 1947) in 1967, when the latter arrived at the Weizmann Institute on a fellowship after finishing his undergraduate studies in physics at King’s College London that same year.³⁶ Warshel and Levitt’s task within Lifson’s group was to write a computer program that would allow calculation of the potential energy of biomolecules viewed as classical spherical atoms bound by harmonic springs whose elastic constants could be evaluated on an empirical ground.³⁷ We now refer to this approach as molecular mechanics (MM), an approximation that allows the BO potential energy of a given molecular system to be expressed as a set of parameterized functional forms that enable the calculation of the structure of the molecular system at its minimum energy.

At the time of Karplus’ visit to the Weizmann Institute, as we pointed out earlier, Lifson’s group was a pioneering reference for the scientific community in the study of potential energy surfaces for small molecules and was a leader in the use of force fields (FF) for complex molecules.³⁸ The group gained great visibility thanks also to the development by Lifson and Warshel of the so-called consistent force field (CFF), a FF optimized by direct comparison with experimental data regarding the geometries and vibrations of the molecules under investigation.³⁹ “In the fall of 1966,” remarked Warshel, “I started my PhD trying to develop what became known as the consistent force field (CFF).”

My general suggested direction was to represent molecules as balls and springs (which became known as molecular mechanics [MM] or a “force field” approach) and to reproduce energies, structures, and perhaps vibrations. This was supposed to be done by a consistent refinement of the MM parameters that will force the calculated and observed properties to be as close as possibly to each other. However, we had no clue how to actually do so.⁴⁰

Yet, the two researchers soon found a way out, or rather, a way into the complexity of intramolecular atomic interactions. Their insight was to introduce in their calculations nonbonded interactions between the atoms of the same molecule—a crucial ingredient that had not been taken into account until then. In particular, they included van der Waals and electrostatic interactions among the partial charges attributed to atoms that are not directly bonded covalently. This was a significant innovation because previous calculations of vibrational molecular properties traditionally considered only the elastic constants governing the oscillations of the atoms around their equilibrium positions to match experimental infrared spectra. The perspective Lifson and Warshel were able to develop was very versatile. In the context of CFF, it is always possible to vary the functional forms to obtain better agreements with the experiments, and the method was immediately applied to organic compounds, thus allowing to expand the applicability of MM (which until then worked for a limited set of small atomic aggregates). Indeed, as we can read in the paper that Lifson and Warshel published in 1969, CFF is “an inductive method for a systematic selection of energy functions of interatomic interactions in large families of molecules.”⁴¹ These developments sparked the curiosity of Karplus who, having masterfully developed techniques for computing the BO surface for ${\rm H}_{3}$, was still essentially unfamiliar with MM and CFF. However, as he explained, “Lifson’s group developed CFF to study simple molecules so as to be able to determine their structure, as well as their vibrational frequencies. They were not interested in biology.”⁴²

As early as December 1969, promising applications of CFF were already beginning to emerge with the publication of a paper by Levitt and Lifson on the refinement of the molecular structure of myoglobin and lysozyme.⁴³ Their study involved CFF to implement a method for refining the Cartesian coordinates of a macromolecule toward its most stable conformation by using an energy minimization procedure. The Weizmann Institute was therefore the right place to be for those interested in these emerging innovations. Concurrently, it was there that Karplus met biochemist Christian Anfinsen (1916–95), who was at the Institute for a visit in 1969.⁴⁴ It was during this meeting that Karplus realized that the protein folding mechanism represented a feasible and compelling intellectual challenge that would allow him to extend to larger and more complex molecules the dynamical approach to scattering reactions he co-developed in the mid-1960s. Indeed, as Karplus recalls:

When I came to Harvard in 1966 and worked on the (H, ${\rm H}_{2}$) reaction, I was not yet thinking of applying the same method to macromolecules. When I went to the Weizmann Institute, I did not know that I was going to work on more complex molecules. It was the lectures of Anfinsen and the CFF potential that Lifson’s group had developed, which gave me the idea that I could apply the method to larger molecules, including biomolecules.⁴⁵

Specifically, Anfinsen, who was shortly to receive the Nobel Prize in Chemistry,

[…] described the experiments that had led to the realization that many proteins can refold in solution, independent of the ribosome and other parts of the cellular environment […]. What most impressed me [Karplus] was Anfinsen’s film showing the folding of a protein with “flickering helices forming and dissolving and coming together to form stable substructures.” The film was a cartoon, but it led to my asking him, […], whether he had thought of taking the ideas in the film and translating them into a quantitative model [italics added]. Anfinsen said that he did not really know how one would do this, but to me it suggested an approach to the mechanism of protein folding, at least for helical proteins such as myoglobin.⁴⁶

It was this idea of a quantitative model, well connected to the article by Hubbard and Wald that appeared in the volume celebrating Pauling’s sixty-fifth birthday, that would begin to shape Karplus’ thinking. As we will discuss in the following section, that essay also played a role during a discussion, which took place at the Massachusetts Institute of Technology (MIT) in 1971, between Karplus and Max Perutz (1914–2002) on the topic of hemoglobin and which helped Karplus engage in the development of an allosteric model to explain its cooperativity. We can thus read Karplus’ early years at Harvard, those between 1966 and 1970, as a period in which he tackled a variety of topics, but without really breaking through as he had done before. It was during this period, in Israel, that he found the spark that would help him restart new impact studies:

During my visit to Lifson’s group, I learned about their work on developing empirical potential energy functions. The novel idea was to use a functional form that could serve not only for calculating vibrational frequencies […] but also for determining that structure. The “consistent force field” (CFF), introduced by Lifson and his coworkers, included nonbonded interaction terms, so that the minimum-energy structure could be found after the energy terms had been appropriately calibrated […]. The possibility of using such energy functions for larger systems struck me as potentially very important for understanding biological macromolecules such as proteins, though I did not begin working on this immediately.⁴⁷

It is useful to clarify, however, what motivated Karplus as he sought to reorient his work toward biology. “It was not that I was dissatisfied with research I was doing during 1966–69 and that I went through a time of creative difficulty,” he explained to us, “but rather that I was looking for a biological problem where what I know how to do would be applicable. That I found in retinal.”⁴⁸ The right path had been found and, as in the case of reactive scattering, he had the insight to find the right collaborators to start moving into a field that was fundamentally new to him. “I went to the Weizmann Institute specifically to switch to biology,” he told us in an e-mail exchange. “Where I did have difficulties, because no one believed it was possible or interesting,” he stressed, “was in extending the methodology used to study the H + ${\rm H}_{2}$ reaction to biomolecules.”⁴⁹

4 Retinal and cooperativity studies in hemoglobin

The only biological subject in which Karplus could consider himself prepared to tackle in 1969 was retinal—a topic he had explored both by working with Hubbard and Wald as a student at Harvard College and by some, alas unsuccessful, attempts at the beginning of his doctoral work at the California Institute of Technology under Delbrück.⁵⁰ He certainly had no institutionalized training in the subject, but, drawn to the topic of vision, he had the opportunity to study retinal, both theoretically and practically.⁵¹ Let it be clear that retinal isomers are linear polyene aldehydes, many of whose physical properties are determined primarily by their π electrons—more delocalized electrons that are often directly implicated in the properties of excited states.⁵² Retinal is a polyene chromophore of visual pigments and is found in both rod and cone cells in the retina of vertebrates, insects, and crustaceans.⁵³ One of its isomers, the one that is of interest to us in this article, is 11-cis-retinal, and it was Wald who understood its function.⁵⁴ Under exposure to visible light, retinal can be found isolated; in the dark, however, it threads its way into the binding site of a protein called opsin, forming rhodopsin, which is a visual pigment. This is the chemical reaction underlying vertebrate vision. In a conversation with us, Karplus remarked that he was very confident about starting a retinal research project and that it was the aforementioned article by Hubbard and Wald led him to the photoisomerization of 11-cis-retinal. The latter, however, was too complex a molecule to be treated with the FCS method that Karplus had co-developed to deal with the study of scattering reactions, but there was hope with CFF semiempirical potentials. As he remarked: “I realized […] that polyenes, such as retinal, were ideal systems for study by the available semiempirical approaches. If any biologically interesting system in which quantum effects are important could be treated adequately at that time, retinal was it.”⁵⁵ Yet, there was a gap to fill and clearly more work to do. The original CFF method, in fact, did not work for chemical reactions because it was not able to take into account the breaking and formation of intermolecular bonds. The task of perfecting the method was challenging, but Karplus had two important allies at his side: Honig, who was already working for him as a postdoc at Harvard, and Warshel, whom Karplus would hire, and who was committed to expanding the CFF method. As the latter noted, in fact:

[…] I started a postdoctoral position at Harvard with Martin Karplus at the beginning of 1970, hoping to make the QM + MM CFF more general. Karplus and his postdoc Barry Honig were at that time making important advances in the study of retinal (the chromophore of the visual pigment), which involves a 12π-electron system. This seemed to be a good rationale to start developing the CFF for π-electron systems. Indeed, I succeeded to connect the molecular orbital (MO) description of atoms with π-electrons with an MM description of σ-bonds with localized electrons, and to consistently refine the corresponding parameters for a unified CFF description. This QM(MO) + MM model included only the bonding between the QM and MM region, and thus ignored all key (e.g. electrostatic) coupling between the MM and QM regions. Nevertheless, the model provided a very powerful and general way to treat large conjugated molecules. During this project, I also figured out how to get the exact analytical forces from the QM treatment, by fixing the molecular orbitals and differentiating only the integrals. As usual, I made this fundamental advance by guessing it, then (as before) I confirmed my idea by using numerical derivatives and then finding the exact mathematical proof. Here again it was shown that the combination of intuition and numerical validation is a powerful tool.⁵⁶

This was an exciting project for Warshel, because the CFF method was inadequate to handle π electrons; it was also an interesting topic for Honig, who had found an appealing ground on which to launch his career under the supervision of a scientist he held in high esteem. And Karplus, for his part, had hired a team that would allow him to land in the much-coveted biology. In other words, the group was well-coordinated, composed of talented researchers and within which everyone had something to gain.

Karplus returned to Harvard in early 1970, and the new team set to work. The first notable results came from the collaboration between Honig and Karplus. As we can read in the introduction of the communication that the two managed to publish in 1971:

Because 11-cis retinal is the chromophore of the visual pigment in vertebrate rods, and its isomerization to the all-trans form […] is known to occur during visual excitation, the properties of the retinal isomers are of considerable interest. Although it has long been assumed that non-bonded repulsions in 11-cis retinal distort the polyene chain, the details of the ground state geometry have not been determined, nor has the experimental spectrum been correlated unequivocally with the excited states.⁵⁷

Honig and Karplus proposed to set up a purely theoretical study, in which they calculated the ground-state potential energy function “with a form which may help to explain some of the unusual spectral characteristics of the retinal isomers when in solution or incorporated into the native visual pigment.”⁵⁸ Using an IBM 7090 they theoretically predicted a possible three-dimensional structure of the visual chromophore. It should be noted that these results appeared, albeit not without difficulty, in Nature, and this was the first time Karplus had published in that journal. This circumstance was undoubtedly a great reward, a motivation that helped fuel the belief that his decision to pursue research in biology was a promising path on which to proceed.⁵⁹ Warshel, for his part, was committed to implementing an optimized and generalized approach of the CFF method to conjugated molecules.⁶⁰ “A detailed interpretation of the electronic transitions and photochemical processes in the retinal isomers requires a knowledge of ground- and excited-state potential surfaces” and the team’s efforts were directed toward that end.⁶¹ In an article published in 1972, Warshel and Karplus provided compelling examples of the applicability of the CFF method to various conjugated molecules, and in 1975, the two published another paper describing a semiclassical approach to photoisomerization, in which classical trajectories were calculated, similar to those used to approach the scattering reaction (H, ${\rm H}_{2}$), in the case of cis–trans isomerization from the triplet state of 2-butene.⁶² The authors, in particular, reported results of “semiclassical trajectory calculations that include non-adiabatic effects and treat the motion in the full multidimensional cartesian space of the molecule undergoing the isomerization reaction.”⁶³ Such retinal trajectories represented, in effect, the first MD simulations of a biomolecule by the Karplus group, although not of a protein^.64 The 1969 classical molecular mechanical CFF model developed by Lifson, Warshel, and Levitt, and which was expanded in 1972 by Warshel and Karplus to incorporate the effect of electronic transitions, was later improved by Warshel and Levitt, culminating in 1976 with a trailblazing article that showed how to combine quantum mechanical (QM) and MM models together to study enzymatic reactions—a computational strategy that is now known as QM/MM simulations.⁶⁵

In the course of his retinal studies, Karplus also undertook another line of research, aligning himself with Attila Szabo (born 1947), a young graduate student originally from Hungary. Emigrating to Canada with his parents in 1957 following the 1956 Hungarian Uprising, Szabo completed his undergraduate studies in chemistry in 1968 at McGill University and then went to Harvard to pursue his doctorate under Karplus. The two produced a simple discrete-state statistical mechanical model of the mechanism of cooperative oxygen binding to hemoglobin, a topic that, by the early 1970s, had been the subject of interest for at least sixty years.⁶⁶ Their goal was not to approach the allostery of the hemoglobin tetrameric protein from an atomic perspective, as had been attempted before them (e.g., Perutz’s stereochemical mechanism), but to represent the nature of the cooperative mechanism through a simple model comprising a finite number of key functional states—a kind of model that would later serve as a blueprint for mapping all-atom MD simulations in terms of dominant metastable states governing biological allostery.⁶⁷ Although Perutz had addressed this topic in qualitative terms, “a quantitative interpretation of the available equilibrium and kinetic data was not attempted.”⁶⁸ More precisely, in fact, in his studies Perutz had already determined the X-ray structure of both oxyhemoglobin and deoxyhemoglobin, and was working on a qualitative explanation of the cooperativity mechanism.⁶⁹ The problem remained open and was very topical.⁷⁰ The underlying assumption was that the hemoglobin molecule exists in two possible quaternary conformations, T (tense) and R (relaxed), that there are two tertiary structures, bound and unbound, for each of the subunits, and that the coupling between the two is ushered in by certain (pH-dependent) salt bridges whose existence is contingent on both the tertiary and quaternary structures of the molecule.⁷¹ The line of research undertaken by Szabo and Karplus aimed to translate Perutz’s mechanism into a statistical mechanical model that could be used to make quantitative predictions to be compared with experiments.⁷² Their first paper appeared in 1972.⁷³

Karplus was then on sabbatical at the Université Paris-Sud (Paris-XI) in Orsay.⁷⁴ “I also began working on hemoglobin,” underlined Karplus in a 2017 interview:

With Attila Szabo, I developed a statistical mechanical model of hemoglobin that transports oxygen in the blood from the lungs to the tissues. I wrote most of the resulting paper in 1972 in Paris at Les Deux Magots, a Left Bank cafe associated with artists and intellectuals, from Hemingway to Sartre. I was on sabbatical leave in Paris, and since 1970 had been planning a permanent life in France.⁷⁵

He then collaborated with the research group of Jeannine Yon-Kahn (1927–2020), an Orsay enzymologist who would become one of France’s leading experts on protein folding and dynamics, developing MD and modeling approaches.⁷⁶ During his sabbatical, Karplus spent a good deal of his time at the Institut de Biologie Physico-Chimique—a research institute, officially separated from Orsay, aimed at promoting research in all areas of biology.⁷⁷ Szabo worked with him at the Institut. In their studies appears a Boltzmann-weighted sum over all key functional states, i.e., the partition function that defines the free energy in statistical mechanics. The model had a definite “Ising model” flavor—a crucial problem in statistical mechanical studies of phase transitions, relying on similar symbolic diagrammatic expansions.⁷⁸ The authors formulated a “generating function to express the relative stabilities of the structures associated with the oxy and deoxy quaternary conformation of the hemoglobin tetramer,” and they introduced “diagrams that display the physical content of the generating function.”⁷⁹

Once Attila Szabo had finished the statistical mechanical model of hemoglobin cooperativity, I realized that his work raised a number of questions that could be explored only with a method for calculating the energy of hemoglobin as a function of the atomic positions. No way of doing such a calculation existed. Bruce Gelin, a new graduate student, had begun theoretical research in my group in 1967.⁸⁰

Relating the effective free energy states from the Szabo-Karplus hemoglobin model to an atomic model of the protein required the ability to determine the energy associated with the atomic coordinates. This was a necessary step, but not a sufficient one, however. Bruce Gelin (born 1946) became the reference figure on whom Karplus could rely on to develop what he had in mind and help him deal with the new challenges that had arisen.⁸¹ Meanwhile, Honig, who had left Harvard in 1970, had completed a postdoctoral fellowship at Columbia University in 1973, where he worked under the supervision of molecular biologist Cyrus Levinthal (1922–90), whose research at that time concerned protein folding.⁸² Szabo received his Ph.D. from Harvard in 1973 and, the following year, he moved to Indiana University Bloomington.⁸³

5 The turning point: BPTI realistic MD simulations

Retinal and hemoglobin studies both raised the need for a treatment of the conformational flexibility of larger systems via some computationally efficient potential energy surface. This circumstance became most strikingly evident with the work on the hinge-bending mode in lysozyme performed by J. Andrew McCammon (born 1947), who was recruited a postdoctoral fellow in the Karplus group in 1976.⁸⁴ The results of his research appeared the same year in Nature, which had led him to produce a harmonic oscillator model for the opening and closing of the active site of lysozyme, and a dynamical model for the overdamped fluctuations of the active site.⁸⁵ The paper was also signed by Gelin, Karplus, and Peter Wolynes (born 1953), a physical chemist who earned his Ph.D. at Harvard and was then a postdoc at MIT with John Deutch (born 1938). “I don’t recall exactly how the work on lysozyme began,” admitted Gelin in a conversation with us:

but my recollection is that I was the one who saw the hinge in a certain drawing and had the idea of imposing an axis around which to open and close the hinge. This was actually a very quick calculation: it involved simply rotating one part of the molecule about the defined axis and doing a single energy calculation at each rotation value. We then refined the process by minimizing the energy at each rotation value (allowing the molecule to ‘relax’ and remove close contacts that gave artificially high energies). The computational work was probably done in just a week or two. The key contribution from Andy [McCammon] was figuring out how this motion would take place within a medium—would it be free, damped, or overdamped? This is where his work with John Deutch complemented the purely computational approach.⁸⁶

McCammon had just finished his Ph.D. thesis at Harvard in 1976 on the subject of hydrodynamics of biopolymers under the guidance of Deutch and Karplus.⁸⁷ He began his doctorate at Harvard in 1969 but had to interrupt it for two years in connection with the Vietnam War.⁸⁸ In 1972 he returned to graduate school, but things went differently than he might have imagined, since, during his doctoral studies, Karplus had the idea of moving to France:

Upon returning to the Ph.D. program in 1972, I started a project on computational study of protein conformations with Martin Karplus. Martin soon decided to relocate to Paris, possibly a permanent move. I approached John Deutch, who was on sabbatical from MIT at Harvard; John kindly agreed to become my day-to-day thesis supervisor, with Martin as cochair of my thesis committee. I proceeded to work on a simple model of phase transitions in lipid bilayers, and more extensively on theoretical hydrodynamics of biopolymers, completing my Ph.D. in 1976. In the process, I had accidentally learned the basics of equilibrium and nonequilibrium statistical mechanics, which proved handy in later work.⁸⁹

Although Karplus’ efforts to open a laboratory in Orsay were partly in vain, his second stint in France gave him the opportunity to think about new ways forward. McCammon was one of the emerging players in efforts toward atomic simulations of proteins, a promising element to put to work with Gelin, who was a driving force in the group.⁹⁰

In his doctoral dissertation, Gelin set out to develop a program that allowed one to work on a given amino acid chain by calculating its potential energy as well as the derivatives of that energy as a function of the positions of the atoms within that particular amino acid sequence.⁹¹ This energy could then be used to find a new stable structure for the polypeptide chain by minimizing its own energy.⁹² Gelin did his thesis work using the Columbia University computer, graduating in 1976 and helping to restructure Warshel and Levitt’s program, which was not a particularly user-friendly code. The computer available at Harvard involved out-of-pocket costs, and Karplus wisely maintained his ties to Columbia after finishing his term there.⁹³ “Earlier Karplus students,” observed Gelin in discussion with us, “worked on 7090 computers at Columbia and Harvard, but most of my work while in the U.S. was done on an IBM 360 Model 91 mainframe at Columbia.”⁹⁴

This was a rare bird: IBM manufactured only a few of them, mainly to hold off competition from Control Data Corporation, whose machines were noted for their scientific calculation capabilities. The 360/91 was sometimes referred to as the Bugatti of computers, because when it was working, it was fantastic—but it was in the shop (being repaired) a lot! Harvard’s Chemistry Department had a remote job entry station just off the lobby of the chemistry building’s Oxford Street front door. There were two keypunches, a card reader, and a line printer, so students could submit jobs and print the results. There must have also been some device to allow one to query the job progress and know when it was done, but I can’t remember what that was.⁹⁵

Gelin achieved substantial results even before he finished his doctorate. Already in 1975, he studies the motion of amino acids sidechains within BPTI, co-publishing with Karplus an article in the Proceedings of the National Academy of Sciences of the United States of America.⁹⁶ This achievement, in which Gelin played the main part, was obtained not only because of Karplus's guidance, but also because of Lifson’s earlier work, as well as the help of Warshel, who made his CFF program available to the Harvard lab. Furthermore, Gelin was able to demonstrate how the weakening effect of heme induced by oxygen binding was transmitted to the interface between hemoglobin subunits: the theoretical analysis that Gelin and Karplus set up provided new insights into the cooperative allosteric mechanism by detailing, at the atomic level, how communication between protein subunits occurs.⁹⁷ These studies related conceptually to the statistical mechanical model of Szabo and Karplus through an explicit atomic representation of the protein. Although all of Gelin’s contributions were purely computational, there were still no MD simulations. Nevertheless, the project of treating an entire protein was beginning to take shape, but it was complex to take the first steps. “I actually considered getting a very large piece of paper and drawing out the whole protein, extended atom by extended atom,” Gelin explained to us:

Fortunately, I thought about the problem a little and realized that I might one day have to define another protein, so I’d better come up something a little more automatic. That led to the idea of the “Residue Topology File” or RTF with the atoms numbered in such a way as to take account of how each amino acid links to the next one. This way, one only needs to specify the amino acid sequence to have the lists of bonds, angles, torsions, and nonbonded interactions enumerated. Then there have to be a couple of custom “patches” applied to take care of things like disulfide cross-links, water molecules, and non-amino acid groups (such as the heme in hemoglobin). The final resulting file I named the “Protein Structure File” or PSF. I was greatly amused years later to see that the RTF and PSF names have been generalized and are now used in building any molecule for computation with CHARMM.⁹⁸

“To implement these ideas,” Gelin, continued:

I created a sequence of small programs (because small programs got higher priority in the batch-processing real-memory computer systems of the day) that passed data from one to the next by files on hard disk storage. To do that I had to learn some aspects of IBM’s Job Control Language (JCL), a barrier that vexed many a student. Using JCL and knowing about various disk file details in the IBM world, you could store your programs as card images in disk files, and not have to read in large decks of cards, with all the risks of card reader errors, decks with cards out of order, and the disastrous dropped deck. The stored card images could be edited and backed up, making development a little simpler.⁹⁹

Together with Gelin, McCammon began working “on a software to simulate the unconstrained dynamics of all of the atoms in a protein—a ‘molecular dynamics’ (MD) simulation.”¹⁰⁰ The protein chosen was BPTI. “We chose this protein because it was small,” remarked Karplus: only 58 residues and only 458 (pseudo) atoms in the extended atom model.¹⁰¹ Moreover, this was one of the few proteins at the time for which a high-resolution crystal structure was available.¹⁰² The calculation started from a potential function created by Gelin, who, by means of empirical energy functions, examined the nature of the potentials that hold BPTI side chains in their equilibrium positions and compared them with those existing for free amino acids.¹⁰³ “Prof. Rahman, who in 1972 had published the first MD simulations of a realistic atomic model of liquid water with Frank Stillinger, came to visit Prof. Karplus and pointed out that since we had a potential function governing the forces on each atom of our protein system, we could have done molecular dynamics,” remarked Gelin.¹⁰⁴

Martin then mentioned this to me and Andy McCammon. In an undergraduate course, I had learned about Rahman’s simulations with van der Waals liquids, so I knew vaguely what we were talking about. While Andy McCammon prepared a simple integrator – not the more sophisticated one Rahman supplied later – I modified one of my programs by removing (or “dummying”) all the calculations except those for bonds. I then defined a very simple molecule: ${\rm H}_{2}$ – just two atoms, one bond, and no other terms to calculate. I remember doing the first run, with an added instruction causing the program to print out the interatomic distance at each step and looking at the printed results (remember: no graphics in those days!). Lo and behold, the ${\rm H}_{2}$ molecule vibrated!¹⁰⁵

At this point, starting with the forces calculated by Gelin, McCammon managed to create a version of the program that allowed him to calculate atomic trajectories within the BPTI molecule.¹⁰⁶ This means that it became possible to use a computer to study atomic movements, simulating the internal dynamics of the biomolecules in question and the conformation the protein could acquire at a given temperature. The goal was to have a more detailed atomic model of BPTI than the simplified coarse-grained model previously utilized by Levitt and Warshel to study protein folding using energy minimization.¹⁰⁷

Among all the challenges the project posed was the more prosaic problem of getting enough computer time to run the simulation itself. “In the mid-1970s,” Karplus emphasizes, “it was difficult to obtain the computer time required to do such a simulation in the United States; the NSF centers did not yet exist. However, CECAM (Centre Européen de Calcul Atomique et Moléculaire [italics added]) in Orsay, France, directed by Carl Moser, a person with an unusual vision for the future of computations in science, had access to a large computer for scientific research.”¹⁰⁸ Karplus refers here to what is known in English as the European Center for Atomic and Molecular Computing, a scientific organization dedicated to the promotion of fundamental research on computational methods applied to condensed matter physics and chemistry.¹⁰⁹ Currently in Lausanne, in the 1970s CECAM was located in Orsay, at the Université Paris-Sud (Paris-XI), where it was founded in 1969 by chemist Carl Moser (1922–2004), who was its first director.¹¹⁰ CECAM was located in the same building as the Centre Inter-Régional de Calcul Électronique (Interregional Centre for Electronic Computing, or CIRCÉ) and could make use of the computational resources provided by the Center.¹¹¹ The year McCammon finished his doctorate, which was 1976, an eight-week extended workshop was organized by biochemist Herman Berendsen (1934–2019) at CECAM, from May 24 to July 17. The topic covered was “Models for protein dynamics” and the idea was to admit only a small number of participants highly specialized in the subject matter:

A number of pioneers in the studies of protein computer simulation were present. Shoshana Wodak, for instance, had completed an algorithm for the analysis of both protein-protein docking and the docking between a protein and other small organic molecules; Michael Levitt was developing a simplified representation of proteins; Peter Rossky and Andrew (Andy) McCammon, of the Karplus group at Harvard, and Wilfred Gunsteren, of the Berendsen group in Groningen, were engaged in specific studies of protein molecular dynamics. This was the right time. Everybody arrived at Orsay with the latest news and a strong desire to obtain new results during the course of the workshop.¹¹²

“Realizing that the workshop was a great opportunity, perhaps the only opportunity, to do the required calculations,” Karplus adds, “Andy McCammon and Bruce Gelin worked very hard to prepare and test a program to do the molecular dynamics simulation of BPTI.”¹¹³ Their job was to create a program combining the protein MM code with a dynamical propagation algorithm, which was working by May 1976. Gelin did not attend the workshop, but the two of them worked together to get the first rudimentary dynamics program up and running, and then it was McCammon who refined it so that he could run the calculations for the BPTI simulation.¹¹⁴ Rahman provided him with the code for the predictor–corrector algorithm to integrate Newton’s equation of motion and what remained for McCammon to do was essentially to use CECAM’s powerful computer to perform all the necessary calculations. The latter was an IBM System/370 Model 168: it was a mainframe and did not offer screens to work with but only a card-reading machine to submit programs and a dot matrix printer available to the public to take the output.¹¹⁵ There was nothing graphic about it other than the command console used by the machine operators (Fig. 1).

Fig. 1

Source: cds.cern.ch/record/750893?1n=en. Courtesy of CERN

An IBM System/370 Model 168. In this photo we see a model in operation at CERN in September 1977.

¹¹⁶

The spirit of the workshop was meant to be a crossroads of ideas aimed at building realistic atomic models of proteins and calculating their dynamics. “The CECAM meeting was very important in my career,” McCammon recalls:

Beyond the MD work, I met Don Ermak, with whom I developed Brownian dynamics simulation methods for diffusional simulations. I also met Herman Berendsen, a wonderful person who had many far-sighted ideas about simulations; Wilfred Van Gunsteren, just starting his postdoc with Berendsen; Michael Levitt; Anees Rahman; and many other stellar scientists.¹¹⁷

The workshop was divided into four main topic areas: stochastic dynamics, accurate dynamics on biomolecular systems, approximate methods on the structure and dynamics of bio-macromolecules, and finally a small working session on potential functions for proteins. As one of the authors recalls: “We went to work every day. There were several rooms with multiple tables, and everyone got a table and worked there when they weren’t discussing with others. Then there was more than one room for group discussions.”

We would study articles, do analytical calculations, have discussions, prepare code in FORTRAN to introduce something into an existing code or to create one from scratch. This was all work that was done at the table, usually in groups of two. For discussions of three or more people, we normally used blackboards. When the FORTRAN to produce a code had been written, one went down to the big CIRCÉ user room and on a machine, one produced the punched cards that then constituted the package/deck code to be submitted to the card reader.¹¹⁸

Moreover, the analysis of programming errors was carried out completely in the self-service room which, in turn, was attached to the machine room where, in addition to the computer, there were tape readers, other printers (now inconceivably bulky), and all the necessary equipment for the operation of the computer.¹¹⁹

The self-service room was always full of people coming and going. This room consisted of several areas, each of which housed a specialized activity such as punch card production, card reading, printing, etc. Participants would send their programs to the computer as we prepared them, and the computer would decide when to process them. Once processed, the results would be sent to the self-service printer and participants would pick them up and review them. The programs could take several hours, and the various scientists would generally send them in at the end of the day and pick them up the next day. The self-service room was always full of people coming and going. This room consisted of several areas, each of which housed a specialized activity such as punch card production, card reading, printing, etc. Everything was organized in detail. In addition, the workshop should also be valued for its social aspect. Those weeks were undoubtedly very productive and we worked hard, but there were also more jovial moments when we went out to eat in small groups at local restaurants and made friendships that would last a lifetime.¹²⁰

Not surprisingly, then, despite the hard work, McCammon found an atmosphere at CECAM that was conducive to his research. The results of the BPTI simulation appeared in 1977 in an article published in Nature by McCammon, Gelin, and Karplus.¹²¹ Based on an MD simulation of BPTI of 9.2 ps, the authors could describe the magnitude, correlations, and decay of fluctuations on the mean structure (Fig. 2).¹²² The simulation showed that the interior of BPTI is flexible and that the motion of the atoms was fluid-like at ordinary temperatures, as the local movements of the atoms have a diffusional character. These results contrast starkly with the rigid view inferred from X-ray structures: “[t]he extent of the protein mobility was, in fact, a great surprise to many crystallographers,” remarks Karplus, “and is one early example of the conceptual insights concerning molecular properties that have been derived from molecular dynamics simulations.”¹²³ This example gives us a way to reflect on the importance of computer simulations in scientific inquiry, which allow for results that would not be obtainable analytically and, as in this case, may be far from what intuition would lead one to imagine. The results obtained by McCammon at the 1976 CECAM workshop were instrumental in initiating MD simulations of the dynamics of complex molecular structures, such as biomolecules. “It is not an exaggeration to say that CECAM has created the molecular dynamics of proteins, and certainly in aqueous environment,” remarked Berendsen.¹²⁴ Also significant are Karplus’ words in describing the computational resources his group found in France: “It would have been difficult,” he noted:

to do the 9.2 ps BPTI simulation in the USA, given the limited computational resources available to academic researchers at the time. Thus, Carl Moser, CECAM, and the Orsay Computer Center played a significant role in making possible the first molecular dynamics simulation of a protein, an essential element in the 2013 Nobel Prize in Chemistry, though not in the Nobel citation.¹²⁵

Although the computer simulation of the BPTI was carried out by McCammon at the workshop, this achievement should be seen as a team effort in which McCammon and Gelin, making good use of Rahman’s suggestions, were guided by their mentor: it was Karplus who identified the research problem and assisted them in developing this work. Gelin uses particularly eloquent words in describing the role Karplus played in guiding his research work and what he assesses to be a remarkable contribution of his doctoral dissertation:

With the wisdom of many years of hindsight, I realize that although we did some good work, probably the most important thing I accomplished was to write it all down, and organize the computer programs clearly, so that following generations of students could build on my work. I probably spent far too much time making my computer programs “pretty,” easy to read and understand, and efficient, but if it helped later students, it’s a good thing. As for a clearly written thesis, much credit goes to Martin Karplus for making me write it and re-write it until enough details were included and explained.¹²⁶

Fig. 2

The peptide backbone and disulfide bonds of BPTI. a) X-ray original structure; b) time-evolved structure after 3.2 ps. In both drawings, alpha carbons are indicated by circles, while sulfides in disulfide bonds are represented by dashed circles. Image taken from McCammon et al. (1977), on p. 586

Also present at the 1976 CECAM workshop was Peter Rossky (born 1950), a graduate student with Karplus who had been working on the development of a sophisticated MBPT for electronic structure before deciding to jump into biomolecular simulations in the course of his doctoral studies. His task was to work with Rahman to produce a simulation of the alanine-dipeptide in a box of water molecules modeled from the ST2 potential developed by Stillinger and Rahman, a work that required going beyond the program written by Gelin and McCammon.¹²⁷ Importantly, simulating a biomolecule in explicit solvent required demanded the imposition of periodic boundary conditions on the system. In contrast, the MD simulations generated by McCammon at the 1976 CECAM workshop were set up considering the protein in a vacuum, and their program did not support periodic boundary conditions. As Rossky puts it:

I felt isolated working as the last Martin Karplus student to do the MBPT. I told him that I would like to do something more in the mainstream of the group for the last project in my thesis. He had been talking to Aneesur Rahman about water and then Martin Karplus suggested that project to me. At that time, Andy McCammon (my office mate!) and Bruce Gelin had already been working for a while on the BPTI MD program. So I started learning the CFF code but didn’t see a water code until I got to CECAM, as I recall. My job at the CECAM meeting was to work with Rahman to merge his ST2 water code with the CFF peptide code as a PBC dynamics code. The minimum image and integrator (a 5th order Gear predictor-corrector) were already in the ST2 code, and I had to be sure everything was properly applied to the peptide. I had no prior knowledge of anything about MD.¹²⁸

Two articles on this topic were published in 1979, although some preliminary results obtained by Rossky and Rahman appeared in the 1976 CECAM workshop report.¹²⁹ For quite some time, their contribution reported the only simulations of a solvated biomolecule surrounded by explicit water molecules and it was not until 1988 that Levitt and Ruth Sharon published a paper detailing a simulation of BPTI surrounded by an explicit aqueous environment of water molecules.¹³⁰ Rahman also collaborated with Jan Hermans (1933–2018) on a study of water dynamics in phase-transition-induced (PTI) single crystals, but although their contribution appears in the CECAM report, it apparently did not lead to any publication.¹³¹ In this regard, it is important to dwell on the technical difficulties that were encountered by these scientists. In many ways, it is not easy for us to understand the challenges these scientists experienced nearly fifty years ago, some of which appear very different from those we typically face today. For instance, Rossky describes his personal challenge to merge the CFF peptide model into the ST2 program of Rahman:

In the early 1970s, people didn’t write codes with the expectation of producing a card deck that would be generally shared. And Anees [Rahman] wrote his own excellent code, but only considering that he should be able to read it. In particular, to save time at the card keypunch, he used a variable naming scheme based on using the fewest keystrokes (x, x1, x2, a, a1, a2, a3, etc…), rather than using names allowing the reader to recognize what variables stored (e.g., xH1, xH2, xO, VxO, FxO, …). And there was no need in such a context for descriptive comment cards.¹³²

Challenges were not uniquely in the code:

Bruce Gelin had already been working for a while on the BPTI MD program, which was not as simple a construction as one might expect (Bruce was an outstanding programmer, and it took a while to get it all debugged). Every parameter in the database for the full polypeptide model was on a separate punch card, and the definition of every term in the potential function for the molecule of interest was also defined on a separate card. It was difficult to avoid all errors.¹³³

Indeed, in the midst of the workshop there were also some mishaps that resulted in amusing anecdotes. For example, McCammon told us that the first MD simulation of a protein he performed was catastrophic: “Martin left it up to me to do what I thought best for the simulation, and my first effort was to examine unfolding. Essentially, the protein blew up within a few picoseconds. Needless to say, this simulation was never published! Feel free to refer to this as a good teaching moment!”¹³⁴ Some recall stories of simulations in which a Lennard–Jones term was missing, causing a remote methyl atom to pass through an adjacent atom almost by magic. Others recall that, after generating one of the first simulations of solvated BPTI in a water box, it was realized that the hydrogen bond between the water and the protein had been inadvertently deactivated. With these stories, we do not mean to pick on anyone in particular, but we do want to shed light on the pioneering nature of this research and the actual difficulty of the work that was being done.

6 Aftermath

As is often the case in scientific research, it is not always possible to know in advance what a particular study may lead to, even after it has been completed; this, incidentally, was the case with the first MD atomic simulations of a protein. In fact, based on our conversations with Karplus, McCammon, and Gelin, the main goal in the years leading up to the publication of the 1977 BPTI article was to make an atomic model of a protein as realistic as possible by relying on MM. At that point, having the ability to calculate atomic forces from that model, the goal was to generate a classical unconstrained trajectory using Newton’s equations of motion. As McCammon notes (Fig. 3), giving us a vivid sense of the enthusiasm of the time:

We did have the sense of doing something new when we tried to take particularly basic equations of motion from physics - very fundamental equations, like Newton’s equation of motion or the diffusion equations - and try to apply them in a detailed way to individual protein molecules, individual nucleic acid molecules. That was something that really had not been done before. It’s something that even at the time we had a sense of bridging in some sort of rather philosophically interesting fashion two really very different areas of science. I mean, here we were, and we had no idea at the time what this might be useful for, if anything. It was purely abstract research at the time. And perhaps a good example of how something that starts out as pure research, not goal oriented in any sense whatsoever, turns out to be tremendously useful and practical for applications later on. We saw this work as an opportunity to use the power of computers, which of course is always a moving marker. We were working with the best machines of the day and with the power of these machines we were able to bridge the basic concepts of physics on the one side with the basic elements of life on the other and to say something about how these molecules actually function. I remember quite vividly when we did the first simulation of the atomic motion in a protein, we did not have nice graphic displays and these things, but we were able to output some of the structures to pen plotters and rather laboriously draw one snapshot of a protein and then a snapshot of what it looked like a little bit later and then a snapshot of what it looked like a little time later. And there was this sense, even at the time, of something truly historic going on, of getting these first glimpses of how an enzyme molecule might undergo internal motions that allow it [to] function as a biological catalyst.¹³⁵

Fig. 3

Photo courtesy of J. Andrew McCammon

Martin Karplus and J. Andrew McCammon (right) in a 1982 photo taken on the occasion of a workshop on biomolecular dynamics organized in Gysinge, Sweden, between October 4 and 7.

In later years, this view was greatly expanded, as MD simulations came to be increasingly perceived as a tool for performing statistical mechanics on biological macromolecular systems. This perspective was already prevalent in fundamental theoretical studies of liquid systems and had begun to consolidate in the 1960s, thanks to the work of Alder, Rahman, Stillinger, Loup Verlet (1931–2019), Bruce Berne (born 1940), David Chandler (1944–2017), Jean-Pierre Hansen (born 1942), John Valleau (1932–2020), and Charles Bennett (born 1943), to name a few prolific authors.¹³⁶ However, statistical physics had not yet been applied in the context of large biological macromolecules. When MD simulations of proteins based on realistic atomic models became feasible, many of the tools and concepts developed in the community performing statistical mechanical studies of liquids could be, and were, readily transferred and implemented in this new field. Some of the most important methods of interest here include the umbrella sampling method used to calculate the relative stability of different conformations, the reactive flux formalism to calculate the transition rate of microscopic processes, as well as the alchemical free energy perturbation (FEP) and thermodynamic integration to calculate the effect of mutations (and chemical changes) on solvation, thermodynamic stability, and ligand binding.¹³⁷ These advanced methods percolated into the field extremely rapidly, often through personal interactions. This was undoubtedly also greatly influenced by the intellectual contributions of Berendsen and Rahman who were both very active participants to the 1976 CECAM. For example, McCammon recalls having the idea of using alchemical FEP to determine changes in drug binding affinity during a seminar by Berendsen describing a thermodynamic study of cavities of different size in water.¹³⁸ This was a direct transfer to the field of biomolecular simulation of computational methods already used in simulation studies of solid and liquids based on a statistical mechanical thermodynamic integration formulation of the chemical potential introduced by John Kirkwood (1907–59) in 1935, as well as a free energy perturbation technique originated by Robert Zwanzig (1928–2014) in 1954.¹³⁹ While it is clear that these contributions provided the fundamental theoretical formulation underlying such calculations, it is equally clear that no one imagined at the time that such formal expressions were meant to be evaluated using MD simulations. As stated by McCammon: “In some sense we were using a set of ideas which were very, very old ideas but just combining them in a completely new fashion.”¹⁴⁰ By the late 1970s, the relationship between statistical mechanics and computer simulations for the study of liquids had matured considerably, and the latter came to be regarded as a legitimate tool for evaluating theoretical expressions formulated through the former.

Interestingly, this profoundly consequential conversion of MD into a computational tool for statistical mechanical studies of biological macromolecules had not really been envisaged before the 1977 BPTI article authored by McCammon, Gelin, and Karplus. As the latter explained: “I was not interested in statistical mechanics per se, but I realized that it was essential for studying biological phenomena at the molecular level.”¹⁴¹ While earlier work, such as the development of the allosteric model of hemoglobin developed with Szabo, had clearly emphasized the centrality of a statistical mechanical framework for understanding protein function, at the time of the 1977 BPTI simulation the MD trajectory of the protein was not perceived as a “tool” for performing statistical mechanics, but primarily as a way to see the “jiggling” of atoms in a realistic model of the protein. All these fundamental ideas simmered just below the surface and, within a few years of 1977, the concept of exploiting MD simulations of proteins to perform statistical mechanics of proteins had been fully assimilated by the biophysical community.¹⁴² In the late 1980s, this view of protein physics matured considerably, and the structural fluctuations that occur within proteins came to be regarded as a key element in understanding their functioning. As Karplus clearly noted at the time: “[i]nstead of viewing protein functions such as enzymatic catalysis only in terms of the structural data provided by high‐resolution X‐ray crystallography, one now recognizes the important role of the internal atomic motions […]” a fact that “opens the way for more sophisticated and accurate interpretations of biomolecular properties.”¹⁴³ In a passage that appeared in 1986 we also read “[t] he molecules essential to life are never at rest; they would be unable to function if they were rigid. The internal motions that underlie their workings are best explored in computer simulations.”¹⁴⁴

Immediately after 1977, user-friendly computer programs and accurate FF development were needed to pursue the goal of rigorously studying protein physics, and both of these tasks were pursued very readily. The suite of programs created by Gelin allowed one to start from the sequence of a protein to create the protein structure file (PSF) that contained the list of all the bonds, angles, and dihedrals within the molecules, as well as all the “types” ascribed to each atom.¹⁴⁵ What Gelin had done represented a great achievement, but the study of any new system still required a lot of careful work with possible human error. To address this issue, the Karplus group’s intensive work led to the creation of the CHARMM program, officially released in 1983, which represented the organized combination of several programs that previously had to be patched together into one integrated MD package.¹⁴⁶ In contrast to this, CHARMM was:

[…] general enough to treat efficiently a range of molecules from simple peptides to the hemoglobin tetramer, and systems from a single molecule to a full crystal composed of a protein and several hundred solvent molecules, while it simultaneously preserves the adaptability and maintainability essential to an academic environment where new projects and application requirements are frequently called for.¹⁴⁷

The basic approach was that of using empirical energy functions for modeling macromolecular systems, namely to “read or model build structures, energy minimize them by first- or second-derivative techniques, perform a normal mode or molecular dynamics simulation, and analyze the structural, equilibrium, and dynamic properties determined in these calculations.”¹⁴⁸ The focus was on “molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments:”¹⁴⁹

For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models.¹⁵⁰

Karplus was initially not so enthusiastic about publishing the force fields used, but in the late 1990s the scientific community lobbied for greater transparency and to make them more publicly known, as demonstrated by Peter Kollman (1944–2001) with the Assisted Model Building with Energy Refinement (AMBER) force field and the widely used simulation program bearing the same name.¹⁵¹ The latter, together with the programs GROMACS, developed by Berendsen’s group at the University of Groningen, and NAMD, developed by Klaus Schulten’s (1947–2016) group at the University of Illinois at Urbana-Champaign, became one of the most widely used engines for protein MD simulations.¹⁵²

The early 2000s saw a remarkable change in the way MD was conceived: during this period, indeed, skilled experimentalists such as John Kuriyan (born 1960) began to rely on MD simulations to advance research on tyrosine kinases, while further MD studies on aquaporin and the potassium channel began to provide significant insight into such systems, also leading to the Nobel Prize in Chemistry 2003 “for discoveries concerning channels in cell membranes.”¹⁵³ Indeed, studies from this period helped to increase the stature and popularity of MD simulations in the biological domain. The all-atom CHARMM force field, encompassing proteins, nucleic acids, and lipid membranes, became the workhorse of most biomolecular simulations from 2000 to 2020, particularly of ion channels and membrane proteins.¹⁵⁴ In 2013, Karplus, Levitt, and Warshel were awarded equal shares of the Nobel Prize in Chemistry “for the development of multiscale models for complex chemical systems.”¹⁵⁵ In describing the achievements of the three scientists, the Nobel committee pointed out how they succeeded in making Newtonian dynamics work side by side with quantum physics by devising methods that combine both.¹⁵⁶ “Previously,” we can read in the Nobel Prize in Chemistry 2013 press release,

chemists had to choose to use either or. The strength of classical physics was that calculations were simple and could be used to model really large molecules. Its weakness, it offered no way to simulate chemical reactions. For that purpose, chemists instead had to use quantum physics. But such calculations required enormous computing power and could therefore only be carried out for small molecules. This year’s Nobel Laureates in chemistry took the best from both worlds and devised methods that use both classical and quantum physics. […] Today the computer is just as important a tool for chemists as the test tube. Simulations are so realistic that they predict the outcome of traditional experiments.¹⁵⁷

Without downplaying the importance of QM/MM, which represents a very rich methodology that is still under intense development today, it is fair to say that many have interpreted the award as a recognition of the importance gained by MD simulations of atomic models of biomolecular systems and the consequent applications that have ensued.

Looking back, we can see that the results Karplus and his group were able to achieve in their approach to biomolecules represented a turning point in a long debate that finds its origin at least with the beginning of statistical mechanics in the late nineteenth century. The transformation affecting biomolecular simulations occurred well after physicists had already focused on nonbiological systems, and it was to the credit of the biological community to have readily incorporated these developments into their own computer simulations. We often speak of “realistic” simulations, as in the case of BPTI; in this context the adjective realistic is often used to mean that the trajectory of atoms follows classical physics, and this was essentially Karplus’ view in the mid-1970s. Eventually, however, the field evolved toward MD calculations that aim to produce realistic results, but in an essentially different way. This is the case, for example, in the aforementioned umbrella sampling or in reactive flux formalism, in which simulations are influenced by algorithmic prescriptions, or in alchemical FEP, in which simulations sample intermediate nonphysical states. More recently, in addition, the statistical mechanical framework of biomolecular computations has been further extended by new machine-learning methodologies from the computer sciences.¹⁵⁸

7 Conclusions

Martin Karplus’ early steps in life science were instrumental in establishing the lines of his academic work in later years. These first studies are situated in a perspective in which computers were conceived as indispensable tools for the joint development of statistical mechanics and biochemistry, making it possible to tackle problems that were unsolvable analytically. As we read extensively in his autobiography, Karplus’ interest in biology dates back long before the 1970s, but until then he had not been able to put it to use. As early as 1971, that is, immediately after his visit to the Weizmann Institute of Science, he had begun publishing studies on 11-cis-retinal conformation and thinking about allostery in hemoglobin. This fact paved the way for the influential studies he would co-develop with Honig and Warshel, on the one hand, and Szabo, on the other, while at the same time giving him greater confidence in the domain of biological studies. Only a few years later, i.e., in 1976, he and his collaborators, notably Gelin and McCammon, were able to extend Rahman and Stillinger’s studies to perform MD simulations of a protein. As we have made clear, Karplus at this point was not yet interested in the statistical mechanics of proteins, but only in their internal motions. However, the 1977 article he co-published is not just a milestone in MD simulations of protein dynamics, but also the occasion that allowed him to start grasping the importance of statistical mechanics in understanding biological processes at the molecular level. This publication, considering also its impact and visibility in Nature, played a key role in setting the stage for the effective merger of computational statistical mechanics and biochemistry.

The events recounted and discussed above show telling characters of the trajectory followed by Karplus, who succeeded in approaching biology, after several failed attempts, only in his mid-40 s. With the equation bearing his name and his studies on reactive diffusion, he could have lived a safe academic life as a respected university professor.¹⁵⁹ Yet he strove to achieve ever-new goals, challenging himself, and sparking threads of research at a time when the biological sciences were experiencing a particularly relevant fervor. “The revitalized idea of molecular dynamics in the late 1970s propagated by Martin Karplus and colleagues at Harvard University sparked a flame of excitement” that continues in full force today fueled by the increasing power of supercomputers.¹⁶⁰ Considering the historical developments reported here, we can only glimpse the complex network of connections, collaborations, and exchanges of ideas that unfolded around Karplus between 1969 and 1977. During this period, but also later with the development of CHARMM, he made use of ingenious collaborators who played a key role in achieving the goals he, as coordinator of a major international research group, set out to achieve. Indeed, one of his great abilities was to put together and manage a laboratory composed of talented, creative, and highly committed young people of the highest technical and intellectual standing. In this sense, Karplus can be credited with bringing together a large group of students and postdocs, giving them a chance to express themselves, identifying research problems, and helping them produce what would become capital of undoubted value to the entire scientific community engaged in MD simulations. In essence, he built a first-rate intellectual and scientific legacy through his teaching, mentoring, and research activities.¹⁶¹ “Karplus was working on a lot of things,” Honig pointed out during a conversation with us, “and he had to lead a lot of people: that’s how it works for a big lab.”¹⁶² As much as this approach has been perpetuated throughout a career spanning 1969 to the present, it has been a trump card that has led Karplus and his team to be trailblazers who have pioneered paths that others had not traveled. The passage, by Ralph Waldo Emerson (1803–82), that Karplus chose as the opening of his Nobel lecture fully captures this spirit: “[d]o not go where the path may lead, go instead where there is no path and leave a trail.”¹⁶³