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

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.


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.
1 Indeed, the re-emergence of classical physics in a world of theoretical chemistry then dominated by quantum mechanics represents a rather peculiar circumstance.
2 In this article we aim to explore specific aspects of this occurrence, pondering how physical through which scientists can achieve the goals of MS. 11 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." 12 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." 13 This approach, devised in 1946 by Stanis law 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  et al. of the Metropolis MC method. 14 This happened at Los Alamos through the use of the MANIAC computer, of whose construction Nicholas Metropolis (1915-99) was the director. 15 In contrast, MD involves "the solution of the classical equations of motion (Newton's equations) for a set of molecules." 16 This task was first accomplished, for a hard-sphere system, by Berni Alder (1925 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  and published them in 1957. 17 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. 18 Rahman was then joined by Frank Stillinger (born 1934) to implement the first MD simulation of liquid water. 19 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. 20 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. 21 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.

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). 22 Until the early 1970s, his research had a predominance of electronic-structure quantum mechanical subjects, e.g., nuclear magnetic resonance (NMR), electron paramagnetic 11 A useful reference for understanding the significance and historical development of MS is Battimelli et al. (2020). See also Macuglia (2022). 12 Allen and Tildesley (2017), on p. 147. 13 Allen and Tildesley (2017), on p. 147. 14 Metropolis et al. (1953). For a useful historical analysis consider Mareschal (2018); Mareschal (2019); Battimelli et al. (2020), particularly pp. 102-110. 15 Here we are referring to MANIAC I, where MANIAC stands for Mathematical Analyzer Numerical Integrator and Automatic Computer. Anderson (1986). 16 Allen and Tildesley (2017), on p. 2. 17 Alder and Wainwright (1957). It should be noted that the paper Alder and Wainwright presented in Brussels eventually appeared in the 1958 symposium proceedings: Alder and Wainwright (1958). Consider also Alder and Wainwright (1959); Alder (2009); Battimelli and Ciccotti (2018); Battimelli et al. (2020). 18 Rahman (1964). See also Battimelli et al. (2020). 19 Rahman and Stillinger (1971); Stillinger and Rahman (1974). See also Allen and Tildesley (2017), on p. 2. 20 Schlick (2010), on p. 426. 21 See, for instance, generalizations of the two-dimensional Ising model in Baxter (1982); Allen and Tildesley (2017), on p. 4, p. 541. 22 There is no space in this article to focus in detail on Karplus' biography. The interested reader may consider Karplus (2020). resonance, along with a major component on reaction kinetics. 23 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. 24 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." 25 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. 26 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 (1914Lifson ( -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. 27 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. 28 "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." 29 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 (1924Hubbard ( -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. 30 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). 31 The former had received his doctorate from the Weizmann Institute in 1968, working under Tel Aviv University physical chemist Joshua Jortner (born 1933). 32 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. 33 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 23 See, for example, Karplus and Godfrey (1966); Shavitt et al. (1968); Godfrey and Karplus (1968); Karplus (1968); Caves and Karplus (1969); Morokuma et al. (1969). 24 Lawler et al. (1967); Purins and Karplus (1969). See Karplus (2020), on p. 109. To account for the many avenues explored by Karplus in his very first years at Harvard, let us not forget the textbook Karplus and Porter (1970). 25 Karplus (2020), p. 109. In this context we must clarify the meaning of the verb "to return": since his childhood, Karplus has shown a strong interest in biology. After some retinal research experiences during college (see footnote 50) and after an unsuccessful attempt to start a doctorate in biochemistry under Max Delbrück , he had to go a long way before "returning" to his central interest. Hereafter we aim to reconstruct and analyze the path he followed. We should also note that beginning in the mid-1970s, Karplus' purely chemical studies became increasingly rare and began to link more closely with biology. Discussing with us, Karplus noted that his last purely chemical study was perhaps Cook and Karplus (1985). 26 It was July 1969 when Karplus arrived at the Weizmann, not the fall of that year as is often stated in the literature. This information was confirmed to us by Karplus and also by Attila Szabo (born 1947), who was his student at the time. 27 Levitt and Lifson (1969); Lifson and Warshel (1969). Lifson was an influential person within the Weizmann Institute, having also served as its scientific director from 1962 to 1966. 28  moved to Harvard to work with Karplus. 34 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. 35 Warshel, for his part, was a graduate student who obtained his Ph.D. in chemical physics under . 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. 36 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. 37 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. 38 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. 39 "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. 40 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 34 We know this information from a conversation with Honig. We also asked him why he had specifically chosen to join the Karplus group out of all the possibilities he had, especially when the latter had not yet done anything in biology or even visited the Weizmann Institute yet. Honig's answer was that he had an inkling that Karplus might be interested in what he was doing: Karplus was already a well-known chemist at the time, not only for quasiclassical scattering studies, but even more so for NMR studies and the equation that bears his name. Honig felt confident that he would have support from the Karplus group. "I gave it a try," he told us, "and things went well." Interestingly, however, there was a dramatic change in topic from the initial intention to focus on photosynthetic processes. In a tone of sympathetic irony, Honig pointed out that back then there was much more flexibility than today in handling research grants. 35 Considering Karplus' autobiography, it might seem that he first read Structural Chemistry and Molecular Biology at the Weizmann Institute. Yet this is not the case. According to Honig's testimony, it was 1968 when Karplus suggested that he read the above article by Hubbard and Wald, and both of them, at the time, were at Harvard. Thus it was that Honig began to devote himself to the study and research of this subject. When Karplus moved to the Weizmann Institute, Honig proceeded very independently, but always under the supervision of Karplus, who identified the problem and pointed him in the right direction. The two talked on the phone about once a month. "Without him," Honig told us referring to Karplus, "I couldn't have done the work." The professional collaboration between Honig and Karplus, in terms of co-published articles, was short-lived, though significant. After his postdoctoral work at Harvard, and after going through other phases of his working career, Honig eventually joined the faculty of Columbia University, where he has currently been teaching and doing research since 1981. A member of the National Academy of Sciences and the American Academy of Arts and Sciences, he is a renowned expert in the field of electrostatic properties of macromolecules and the head of a flagship laboratory in molecular biophysics. 36 In 1968 Levitt had begun his doctorate in computational biology at the Medical Research Council (MRC) Laboratory of Molecular Biology (LMB) in Cambridge, completing it in 1971. He continued as a member of the LMB staff until September 1972 and then obtained a two-year postdoctoral position at the Weizmann Institute. He then returned to Cambridge. From 1977 to 1979 he worked at the Salk Institute for Biological Studies and, in the same year, he moved to Israel. More information about Levitt's biography can be found here at www.nobelprize.org/prizes/chemistry/2013/levitt/biographical/. We invite interested readers to also read Levitt (2001). Regarding Warshel's biography, please refer to www.nobelprize.org/prizes/ chemistry/2013/warshel/facts/. 37 See also Miller (2013), a short, yet effective, article that lays out chronology to explain the rationale for the 2013 Nobel Prize in Chemistry. 38 A force field is typically a potential function composed of simple analytical functions together with a set of empirically optimized parameters that is designed to mimic the quantum mechanical BO potential energy of an atomic or molecular system. In 1969, Lifson's was not the only prominent group in this type of research: also worth mentioning are those of Norman Allinger (1928) and Harold Scheraga (1921. , on pp. 9992-9993. 39 Lifson and Warshel (1969). See also Levitt and Lifson (1969). 40 Warshel (2014), on p. 10021. 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." 41 These developments sparked the curiosity of Karplus who, having masterfully developed techniques for computing the BO surface for 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." 42 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. 43 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. 44 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, 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. 45 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. 46 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 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 41 Lifson and Warshel (1969), on p. 5116. 42 Martin Karplus in conversation with the authors. 43 Levitt and Lifson (1969). 44 Karplus (2020), on p. 123. 45 "I heard Anfinsen's lecture in the summer of 1969 when I was in Lifson's group," remarked Karplus in conversation with the authors. 46 Karplus (2020), on p. 123; Anfinsen (1973).
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. 47 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." 48 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 + H 2 reaction to biomolecules." 49

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. 50 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. 51 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. 52 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. 53 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. 54 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." 55 Yet, there was 47 Karplus (2020), on p. 124. It should be noted that at the Weizmann Institute Karplus also devoted himself to proofreading the manuscript of the textbook Karplus and Porter (1970). 48 Martin Karplus in conversation with the authors. 49 Martin Karplus in conversation with the authors. We would like to point out that as early as 1971 Karplus published papers related to retinal and hemoglobin. In the next section, we will comment on the retinal studies that Karplus undertook. See ; ; Shulman et al. (1971). 50 During his second year at Harvard College, where he enrolled in 1947, Karplus studied under Wald and Hubbard. The encounter with them provided him with a clear understanding that the phenomena of life could be approached at the molecular level; indeed, the young student took a course from Wald entitled "Molecular Basis of Life." As part of the course material, there was a section on vision. Equally interesting was Karplus' collaboration with Hubbard, with whom he had the opportunity to work on the visual chromophore. "My undergraduate research with Wald and Hubbard," stressed Karplus in conversation with us, "was essentially experimental and involved isolating components of the visual pigments." 51 Considering the people with whom the young Karplus had interacted, his retinal training, while not part of a biology degree program, nevertheless constituted a first-rate preparation: in 1967, Wald shared the Nobel Prize in Physiology or Medicine equally with physiologists Ragnar Granit (1900-91) and Haldan Keffer Hartline (1903-83) "for their discoveries concerning the primary physiological and chemical visual processes in the eye." See www.nobelprize.org/prizes/medicine/ 1967/summary/. Hubbard, on his side, about whom Karplus reserves words of praise in his autobiography, was a researcher of the highest value, although she obtained an academic position much later than she would have deserved. Karplus (2020), on p. 67. 52 Retinal is a conjugated molecule, i.e., "a group or chain of atoms bearing valence electrons that are not engaged in singlebond formation and that modify the behaviour of each other." See "Conjugated system," in Encyclopaedia Britannica, online resource, www.britannica.com/science/conjugated-system; Miller (2013), on p. 14; Honig et al. (1975), on p. 92. 53 Gilardi et al. (1971), on p. 187. 54 We should point out that Wald never did any structural work on the 11-cis-retinal. Prior to Wald's discovery, Pauling had evaluated 11-cis-retinal as an unlikely candidate for being the visual chromophore. See Honig and Ebrey (1974), on p. 154. 55 Karplus (2020), on p. 111. 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. 56 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. 57 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." 58 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. 59 Warshel, for his part, was committed to implementing an optimized and generalized approach of the CFF method to conjugated molecules. 60 "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. 61 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, H 2 ), in the case of cis-trans isomerization from the triplet state of 2-butene. 62 The authors, in particular, reported results of "semiclassical trajectory calculations that include nonadiabatic effects and treat the motion in the full multidimensional cartesian space of the molecule undergoing the isomerization reaction." 63 Such retinal trajectories represented, in effect, the first MD simulations of a biomolecule 56 Warshel (2014), on p. 10022. 57 , on p. 558. 58 , on p. 558. 59 . Consider Karplus (2020), on p. 112; to this regard see also Honig et al. (1975), on p. 97. 60 This approach was based on an innovative separation of σ and π electrons, representing the former through empirical potential functions and the latter by means of a semiempirical model, of the Pariser-Parr-Pople type, corrected for nearestneighbor orbital overlap. Warshel and Karplus (1972), on p. 5612, p. 5613. 61 Honig et al. (1975), on p. 95. 62 Warshel and Karplus (1972), on pp. 5622-5625; Warshel and Karplus (1975). 63 Warshel and Karplus (1975), on p. 11. 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. 65 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. 66 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. 67 Although Perutz had addressed this topic in qualitative terms, "a quantitative interpretation of the available equilibrium and kinetic data was not attempted." 68 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. 69 The problem remained open and was very topical. 70 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. 71 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. 72 Their first paper appeared in 1972. 73 64 Retinal studies are not the focus of our paper, and we refer the interested reader to Karplus' autobiography and the scientific articles he co-published on the subject. 65 Lifson and Warshel (1969); Levitt and Lifson (1969); Warshel and Karplus (1972), on pp. 5622-5625; Warshel and Karplus (1975); Warshel and Levitt (1976). 66 As concisely and effectively stated in the Collins Online Dictionary, allostery is "the condition of a protein (such as an enzyme) in which the structure and activity of the enzyme are modified by the binding of a metabolic molecule at a site other than the chemically active one." That said, which is essentially correct and intuitive, in an e-mail exchange with one of us (Benoît Roux), Jean-Pierre Changeux (born 1936) explained in his own words that the generally accepted meaning of allostery is to qualify "indirect interactions between topographically distinct sites mediated by a conformation change." For more details see also Changeux (1961Changeux ( , 1964; Monod et al. (1965); Changeux (2013). In relation to the studies of Szabo and Karplus, consider Szabo and Karplus (1972b). See also Szabo and Karplus (1972a). The latter, though, is an abridged version of the former, and was submitted about a month earlier than the second publication. 67 Szabo and Karplus (1972b), on pp. 163-164. Ushnish and Strodel (2018). 68 Szabo and Karplus (1972a, b), on p. 164. Perutz (1970a, b). 69 Perutz (1970a); Morimoto et al. (1971); Greer and Perutz (1971). 70 As briefly alluded, Szabo and Karplus' research work took shape after a discussion that Karplus had with Perutz at MIT in 1971. The latter was invited by biologist and biophysicist Alexander Rich  to give two lectures at MIT, during which he spoke about the structure of oxidized and deoxidized hemoglobin and proposed a stereochemical mechanism to explain cooperativity in such a protein. During their discussion, which followed the second lecture, and which was held privately in Rich's office, the focus was on the experimental results that Perutz had achieved. On that occasion, Karplus suggested to Perutz, contrary to what the latter had done, that the problem of cooperativity in hemoglobin be approached in quantitative, thermodynamic terms. One understands again the importance of the words that Wald and Hubbard wrote in the mentioned publication in Pauling's honor. It was after that meeting with Perutz that Karplus began to take his next step into biology: "[h]aving been taught by Pauling that until one expressed an idea in quantitative terms, it was not possible to test one's results, I went away from our meeting thinking about the best way to proceed." See Karplus (2020), on p. 115; Perutz (1970a); Morimoto et al. (1971). In an e-mail exchange with the authors, Szabo recalls that Karplus asked him to attend Perutz's lectures and that, shortly afterward, the two set to work on building a quantitative model for hemoglobin. At that time, Szabo had been collaborating with Karplus for a little over a year and had also worked on pseudocontact shifts in low spin iron-porphyrin complexes without, however, producing a publication about it. See also Szabo (2008), on p. 5883. 71 Karplus (2020), on p. 115. 72 Note that, prior to publishing with Szabo, Karplus had already worked on hemoglobin, co-producing a quantum-mechanical model "to interpret the shifts measured by nuclear magnetic resonance in a variety of low spin (s = 1/2) ferric cyanide heme and heme protein complexes." See Shulman et al. (1971), on. p 93. 73 Szabo and Karplus (1972a).
Karplus was then on sabbatical at the Université Paris-Sud (Paris-XI) in Orsay. 74 "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. 75 He then collaborated with the research group of Jeannine Yon-Kahn (1927, an Orsay enzymologist who would become one of France's leading experts on protein folding and dynamics, developing MD and modeling approaches. 76 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. 77 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. 78 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. The two wrote and submitted the aforementioned article before leaving for France. 75 See news.harvard.edu/gazette/story/2017/04/harvards-martin-karplus-looks-back-on-path-to-nobel-prize/. 76 Yon-Kahn was a female scholar in a universe dominated essentially by males A succinct biography of her, although in French, can be read by clicking on the following link: www.insb.cnrs.fr/fr/cnrsinfo/hommage-jeannine-yon-kahn. Karplus spent two sabbatical periods in France, the first at Paris-XI between 1972 and 1973 and the second, at Paris-VII, between 1974 and 1976, serving as professeur associé in the academic year 1975-1976. 77 The Institut de Biologie Physico-Chimique, founded in 1927 and inaugurated in 1930, is of particular interest if for no other reason than the substantial funds on which it could rely. Two important names are behind its history: Edmond James de Rothschild (1845Rothschild ( -1934, philanthropist and creator of the Rothschild Foundation, and Jean Baptiste Perrin , recipient of the Nobel Prize in Physics 1926 and first president of the Institut. For further information regarding the Institut please consider: www.ibpc.fr/en/. 78 Szabo and Karplus (1972a). Consider also (Szabo and Karplus 1973 Karplus lab in mid-1969, only to discontinue research in August 1969 to begin military service. Prior to this interruption, he had worked with members of the Karplus group on the application of the random phase approximation to two-electron systems. Upon his return to Harvard in September 1971, he decided to devote himself to the study of conformational problems in small molecules, and later in proteins, partly because during his military service he had been assigned to a unit dealing with illegal drug use. As Gelin told us, "I ended up (by some sheer chance and good luck) in the U.S. Army Criminal Investigation Division's laboratory at Fort Gordon, Georgia, where I worked as a forensic chemist evaluating evidence seized by investigators. The most common findings were marijuana, heroin, prescription drugs, and some hallucinogens such as LSD. The work did create an interest in structure-function relationships." at that time concerned protein folding. 82 Szabo received his Ph.D. from Harvard in 1973 and, the following year, he moved to Indiana University Bloomington. 83

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. 84 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. 85 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. We report for completeness a work on retinal that Warshel produced upon his return to Israel, but which we will not address in our study: Warshel (1976). 83 At the University of Indiana, Szabo got interested in NMR studies of the internal dynamics of biomolecules. See Wittebort and Szabo (1978). Subsequently, with his graduate student Giovanni Lipari (1953, he developed the influential modelfree approach to the interpretation of NMR relaxation. See Lipari and Szabo (1980); ; . It should also be noted-and this will be useful for the following-that Lipari co-published a paper that is considered to be the first comparison of MD with experiments in the case of BPTI. See . 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. 91 This energy could then be used to find a new stable structure for the polypeptide chain by minimizing its own energy. 92 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. 93 "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." 94 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. 95 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. 96 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. 97 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 Footnote 90 continued to Paris-VII. Additional people who joined us there included Henry Suzukawa, David Munch, and Vincent B.H. Hyunh. There may also have been a French student named Daniel Ballou." 91 In an e-mail exchange with us, Gelin explained that: "Spectroscopic data-vibrational spectra-were used to set the energy parameters of the functions that expressed the resistance of bonds, angles, torsions, and other internal coordinates to deformation from their equilibrium values." Gelin's program, which was based on the CFF program, ended up constituting the first version of the CHARMM program (which we will introduce in the next section). As Karplus remarked in conversation with us: "Having the CFF program, brought to Harvard by Warshel permitted Bruce Gelin to develop the first version of the CHARMM program and apply it to a number of problems." 92 Gelin (1976). 93 See Macuglia et al. (2021), on pp. 3-4. 94 Bruce Gelin in conversation with the authors. "In France," he emphasized recalling his visit to Orsay, "I worked on the IBM 370/165 at CIRCÉ, as you note below. It had a very good self-service setup, with lots of keypunch machines, and again a self-service card reader and a fast line printer." 95 Bruce Gelin in conversation with the authors. "At Harvard," Gelin also stressed, "we had the remote job entry station in the chemistry department (the one I described with self-service card reader and line printer). Also in that space just off the lobby of the chemistry department was a PDP-11/45 (Digital Equipment Corporation) but that was not used for molecular modeling; others used it for computationally smaller tasks. Harvard also had its own computer center a few city blocks from the chemistry department; it had some model of IBM mainframe-that was more expensive to use, so I didn't use it much." 96 Gelin and Karplus (1975). 97 Gelin and Karplus (1977); Karplus (2020), on p. 125. 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. 98 "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. 99 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." 100 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. 101 Moreover, this was one of the few proteins at the time for which a high-resolution crystal structure was available. 102 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. 103 "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. 104 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: 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 H 2 molecule vibrated! 105 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. 106 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. 107 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  Gelin and Karplus (1975). 104 Bruce Gelin in conversation with the authors. The meeting between Rahman and Karplus took place at Harvard, as Gelin explained in a discussion with the authors. "This happened after we had all returned from France," he pointed out. "It was then that Martin gave Andy and me the task of adding a dynamics simulation feature to the existing force-field calculations." 105 Bruce Gelin in conversation with the authors. 106 McCammon (2016), on p. 8058. 107 Levitt and Warshel (1975). 108 Karplus (2014), on p. 9997. applied to condensed matter physics and chemistry. 109 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 , who was its first director. 110 CECAM was located in the same building as the Centre Inter-Régional de Calcuĺ Electronique (Interregional Centre for Electronic Computing, or CIRCÉ) and could make use of the computational resources provided by the Center. 111 The year McCammon finished his doctorate, which was 1976, an eight-week extended workshop was organized by biochemist Herman Berendsen (1934Berendsen ( -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. 112 "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." 113 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. 114 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. 115 There was nothing graphic about it other than the command console used by the machine operators (Fig. 1). 116 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. 117 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 109 There is no official English translation of the acronym CECAM. Although it often reads European Center for Atomic and Molecular Calculations, we believe the word computing is more appropriate because it highlights the close connection between this institution and its mission of using, developing and promoting computational methods. See www.cecam.org.  (born 1943) and her work on protein-protein docking, consider the website of her lab at the Centre for Computational Biology at The Hospital for Sick Children, Toronto, Canada: wodaklab.org/ws. We will have a chance to talk more specifically about Peter Rossky (born 1950) later in this article. In relation to Wilfred Gunsteren (born 1947), we recommend reading Hünenberger et al. (2012). 113 Karplus (2014), on p. 9997. 114 Note that Warshel did not attend the workshop either. Levitt, on the other hand, joined it only in part and Karplus was a visitor. "I participated in the workshop only in part," Karplus explained to us, "because I felt that the important aspect was Andy [McCammon] doing the BPTI simulation and reporting on it." This reinforces the idea, already deeply rooted in us, of a mentor who was able to follow his students without overpowering them, but providing them with research ideas, helping them, and giving them the space they needed to become independent scholars. It should also be mentioned, as stated in footnote 76 , that Karplus was also professeur associé at Paris-VII in the 1975-76 academic year, a status that occupied him and left him less time to participate in the workshop full time. 115 See www.ibm.com/ibm/history/exhibits/mainframe/mainframe_PP3168.html. 116 The reader interested in the technical details of this specific model of IBM computer can refer to Pugh et al. (1991). Consider also Germain (1967). 117 McCammon (2016), on p. 8058. See also Ermak (1975). 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. 118 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. 119 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. 120 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. 121 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). 122 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." 123 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. 124 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. 125 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. 126 Footnote 122 continued he probably used programs that I wrote. Andy and I had a good collaboration as we both worked on getting the molecular dynamics to run." 123 Karplus (2014), on p. 9997, see also p. 9993. 124 Berendsen (1990), on p. 10. 125 Karplus (2016), on p. 3. This piece by Karplus is found in an unpublished report that was compiled in 2016 on the occasion of the CECAM meeting "Models for Protein Dynamics, 1976Dynamics, -2016." This report is untitled and, in this note, we have decided for convenience to title it Models for Protein Dynamics, 1976Dynamics, -2016 , which, however, is the title of the CECAM meeting organized to celebrate the 40th anniversary of the 1976 workshop. Please contact CECAM to read this report. See also Battimelli et al. (2020), on p. 103, pp. 108-109. 126 Bruce Gelin in conversation with the authors.
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. 127 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. 128 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. 129 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. 130 Rahman also collaborated with Jan Hermans (1933Hermans ( -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. 131 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. 132 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. 133 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! 127 Stillinger and Rahman (1974). 128 Peter Rossky in conversation with the authors. 129 Rossky and Rahman (1976); ; . The 1976 CECAM workshop report was never officially published; it was recently archived by CECAM at our suggestion and is currently accessible by clicking on this link: heyzine.com/flip-book/8f2654df2d.html. 130 Levitt and Sharon (1988). 131 Hermans and Rahman (1976). Judging from his involvement, it is clear that Rahman played a decisive role in the 1976 CECAM workshop, pushing to combine advances in MD simulations of liquids with modeling of biomolecules. 132 Peter Rossky in conversation with the authors. 133 Peter Rossky in conversation with the authors. Feel free to refer to this as a good teaching moment!" 134 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.

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. 135 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, and Charles Bennett (born 1943), to name a few prolific authors. 136 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. 137 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. 138 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  in 1935, as well as a free energy perturbation technique originated by Robert Zwanzig (1928) 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." 140 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." 141 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. 142 In the late 1980s, this view of protein physics matured considerably, and the structural fluctuations that occur 136 For a comprehensive collection of landmark articles in the field, the reader is invited to consider Ciccotti et al. (1987). 137 Umbrella sampling is a technique used to improve the sampling of a system via computer simulations by introducing a biasing potential whose effect is then removed in final analysis. On the other hand, the reactive flux formalism is a technique employed to calculate the transition rate in a system via computer simulations by initiating trajectory at the top of the activation free energy barrier. Finally, alchemical FEP and thermodynamic integration are techniques used to calculate free energy difference between two states of a system via computer simulations by progressively transforming the potential energy from one state to the other. Regarding the umbrella sampling method consider Bennett (1976); Torrie and Valleau (1977). In relation to the reactive flux formalism, see Bennett (1977); Chandler (1978); McCammon and Karplus (1979); Northrup et al. (1982). With respect to FEP and thermodynamic integration consider, for example, Postma et al. (1982); Tembe and McCammon (1984); Jorgensen and Ravimohan (1985); Lybrand et al. (1986); Bash et al. (1987a, b); Gao et al. (1989). 138 Allison (1995). 139 Ciccotti et al. (1987), on pp. 80-159;Kirkwood (1935); Zwanzig (1954). 140 Allison (1995). 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." 143 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." 144 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. 145 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. 146 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. 147 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." 148 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:" 149 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. 150 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 (1944Kollman ( -2001 with the Assisted Model Building with Energy Refinement (AMBER) force field and the widely used simulation program bearing the same name. 151 The latter, together with the programs GROMACS, developed by Berendsen's group at the University of Groningen, and NAMD, developed by Klaus Schulten's (1947Schulten's ( -2016 group at the University of Illinois at Urbana-Champaign, became one of the most widely used engines for protein MD simulations. 152 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." 153 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. 154 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." 155 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. 156 "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. 157 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. 158

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. 159 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. 160 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 and1977. 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. 161 "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." 162 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." 163