Victoria Madame

Garcia Sakai

PETERS Judith Madame

Natura in minima maxima.

"Nature is the greatest in the smallest things" -

Résumé en franc ¸ais

Le travail présenté dans cette thèse vise à apporter des informations microscopiques sur la stabilité thermique des cellules bactériennes (dénaturation et croissance cellulaire).

La thèse est divisée en six chapitres, à partir du premier qui contient une présentation des effets de la température sur la matière vivante, introduit le thème de la réponse thermique des bactéries à haute température.

Le deuxième chapitre est consacré au cadre théorique de la diffusion neutronique en physique de la matière condensée, dont le formalisme est introduit conduisant aux principales équations pour les observables pertinentes. La diffusion élastique (EINS), quasi-élastique (QENS) et inélastique des neutrons est également définie. Dans le domaine de l'EINS, l'approximation gaussienne et le comportement non gaussien dû à l'anharmonicité et à l'hétérogénéité dynamique sont discutés, et différents modèles d'hétérogénéité sont considérés. Dans le domaine de QENS, il est dérivé de la manière de décomposer le spectre pour fournir des informations sur les mouvements des atomes à différentes échelles de longueur et de temps : rigide (global), diffusif confiné (comme par exemple, les mouvements d'hydrogène des groupes méthyle "diffusifs en cage") et vibrations interne.

Les méthodes de calcul utilisées dans la thèse sont décrites dans le troisième chapitre. Il contient l'explication des techniques standard de dynamique moléculaire, en mettant l'accent sur les simulations de protéines. Deux scénarios sont décrits : des simulations atomistiques complètes et des simulations à gros grains, y compris l'approche fournie par la méthode Lattice Boltzmann.

Le chapitre 4 présente le résultat principal de ce travail : grâce à l'utilisation de la diffusion quasi-élastique des neutrons combinée à des simulations de dynamique moléculaire, il a été montré qu'il y a un fort ralentissement de la diffusion globale dans le cytoplasme à partir juste en dessous de la température de mort cellulaire soutenant l'idée du déploiement des protéines dans le cadre du processus de dénaturation irréversible. Étonnamment, la fraction de protéines dépliées n'est que d'environ 5%, mais c'est l'effet que cela provoque sur l'ensemble du protéome et son comportement, qui affecte la viabilité cellulaire. Aucune dénaturation catastrophique du protéome ne se produit à la mort cellulaire.

Les derniers chapitres ont été consacrés à deux travaux complémentaires traitant tous deux des effets de la complexation protéine-ligand sur la dynamique des protéines, dans lesquels techniques de diffusion neutronique et simulations de dynamique moléculaire sont toujours couplées.

Cells and bacteria

In biology, an individual cell is the minimal self-reproducing unit of living matter. To distinguishe living from non-living systems, we may consider that all living systems metabolize (i.e. take in resources), grow and replicate. However, there are systems that are clearly not alive who fit this definition -e.g. candle flames can take in fuel, oxidize it, grow bigger fires and light new fires; oil droplets can grow and duplicate, and similar processes can occur also in self-replicating computer codes or in human institutions. To exclude these possibilities, here we follow the definition of Dill and Agozzino [START_REF] Dill | Driving forces in the origins of life[END_REF], considering a system alive when it:

• exchanges energy and matter with the environment;

• grows and replicates independently;

• has ancestry and offspring variation;

• is formed by molecules.

The unity and diversity of organisms become apparent at the cellular level. The smallest organisms consist of single cells and are microscopic. Larger, multicellular organisms contain many different types of cells, which vary in size, shape, and specialized function. The cellular volume can vary over four orders of magnitude, from ⇡ 1.3µm 3 in the case of the small Escherichia coli bacterium, up to ⇡ 10 4 µm 3 for the huge mammalian cell [START_REF] Milo | Cell biology by the numbers[END_REF]. Despite the impressive diversity of living systems, all cells of the simplest and most complex organisms share certain fundamental properties, which can be observed at the biochemical level. In particular, cells of all kinds share a few of structural features [START_REF] Nelson | Lehninger Principles of Biochemistry[END_REF] -see Fig. 1.1. Their boundaries are defined by plasma membranes, separating cells' contents from the surroundings. Plasma membranes are made of lipid and protein molecules, creating a fine hydrophobic barrier around the cell that blocks the free passage of inorganic ions and majority of charged and polar compounds. Transport proteins located in the membrane control the passage of specific ions and molecules -receptor proteins transmit signals into the cell, and membrane enzymes participate in some reaction pathways. Moreover, since the individual proteins and lipids that form the plasma membrane are not covalently bound, the entire structure of the membrane is considerably flexible, allowing changes of the size and shape of the cell. As a cell grows, newly made lipid and protein molecules are inserted into its plasma membrane; cell division produces two cells, each with its own membrane. This growth and cell division occurs without loss of membrane integrity [START_REF] Nelson | Lehninger Principles of Biochemistry[END_REF].

In general, cell reproduction involves the transmission of genetic information to their offspring. Each cell stores this information in the same chemical form, i.e. as double-stranded DNA molecules. For the replication, the cell copies its information by splitting up the paired DNA strands and using each strand as a template for polymerization to build a new DNA strand with a complementary sequence of nucleotides. The same strategy of templated polymerization is used by the cells for the synthesis of proteins. This is achieved through the transcription of portions of the genetic information from DNA into the closely related RNA molecules, which in turn guide the protein production by a more complex mechanism of translation that takes place in the ribosomes [START_REF] Alberts | Molecular biology of the cell[END_REF]. These processes of replication and storage of the genetic information, with the associated proteins, take place in a specific region inside the bacteria. For bacteria and archaea, this region is called nucleoid and is quite irregular, meanwhile, in the case of eukaryotic cells, it is enclosed within a double membrane and is called nucleus. The remaining material enclosed between the plasma membrane and the nucleous (or nucleoid) formes the so called cytoplasm. The liquid phase of the cytoplasm in an intact cell is the cytosol, which does not contain any part of the cytoplasm that is inside the organelles. The set of all proteins in the cell is called proteome.

The effects of temperature on bacteria

Temperature is one of the key environmental factors in microbial life, and it can be used to classify different bacteria depending on the temperature of optimal growth (T OG ) -i.e. the thermal condition at which the bacteria live and thrive [START_REF] Kumar Chandrayan | Hyperthermophile protein behavior: partially-structured conformations of Pyrococcus furiosus rubredoxin monomers generated through forced cold-denaturation and refolding[END_REF][START_REF] Galtier | Relationships between genomic G+ C content, RNA secondary structures, and optimal growth temperature in prokaryotes[END_REF][START_REF] Musto | Genomic GC level, optimal growth temperature, and genome size in prokaryotes[END_REF][START_REF] Nakashima | Compositional changes in RNA, DNA and proteins for bacterial adaptation to higher and lower temperatures[END_REF][START_REF] Ohtani | An extreme thermophile, Thermus thermophilus, is a polyploid bacterium[END_REF], see Table 1.1. Generally, extreme temperatures are problematic for all living systems. In the case of bacteria, low temperatures can induce a decrease in enzymatic activity in the cell, and the velocity of various biochemical reactions, cell metabolism, and the fluidity of biomass membranes is reduced [START_REF] Hébraud | Cold shock response and low temperature adaptation in psychrotrophic bacteria[END_REF][START_REF] Shivaji | How do bacteria sense and respond to low temperature?[END_REF]. If the temperature is further lowered below the freezing point, the water inside the cell condenses forming ice crystals, causing irreversible mechanical damage to the biomass membrane [START_REF] Sandra | Selection of low-temperature resistance in bacteria and potential applications[END_REF]. On the other hand, high temperatures provoke the irreversible denaturation of proteins, affecting the normal cell physiological activities, such as damaging the enzymes involved in the tricarboxylic acid cycle [START_REF] Ventura | How high G+ C Gram-positive bacteria and in particular bifidobacteria cope with heat stress: protein players and regulators[END_REF]. Besides, high temperatures can lead to a loss of integrity for plasma membrane and to the damage of nucleic acid molecules [START_REF] Russel | Lethal Effects of Heat on Bacterial Physiology and Structure[END_REF][START_REF] Teixeira | Identification of sites of injury in Lactobacillus bulgaricus during heat stress[END_REF].

In this thesis we will focus on the problem of the thermal response of bacteria at high temperature.

Thermal inactivation of bacteria

Several factors influence bacterial inactivation due to high temperatures such as the growth conditions and composition, and pH of the environment surrounding the cells [START_REF] Stringer | Thermal inactivation of Escherichia coli O157: H7[END_REF]. Every component of the bacteria (membranes, proteins, DNA, RNA) will be affected to some degree by temperature increase. Thus, it is not an easy task to ascribe a lethal effect to a single alteration within an organism. Nevertheless, some changes are more pronounced than others, some types of damage may be repairable and all depend upon the exact value of the temperature applied.

• Damage to the membrane. In the case of Gram-negative bacteria (such as the E. coli), the outer membrane is quite sensitive to high temperature. In particular, damage to the membrane already occurs when cells are subjected to moderate heat shocks [START_REF] Beverley | Outer-membrane damage in sublethally heated Escherichia coli K-12[END_REF]. Morphological and structural changes altering the permeability of the barrier have been observed [START_REF] Katsui | Heat-induced blebbing and vesiculation of the outer membrane of Escherichia coli[END_REF] -see Figure 1.2 and 1.3. This effect is enhanced in the presence of Tris buffer [START_REF] Tsuchido | Destruction of the outer membrane permeability barrier of Escherichia coli by heat treatment[END_REF][START_REF] Tsuchido | Interaction of the fluorescent dye lN-phenylnaphthylamine with Escherichia coli cells during heat stress and recovery from heat stress[END_REF] However, the cytoplasmic (inner) membrane, that is responsible for controlling the flux of molecules entering and leaving the cell interior, is more resistant to heat shocks. Damage to this membrane has crucial effects on bacteria, causing leaking of intracellular material [START_REF] Russell | Potential sites of damage in micro-organisms exposed to chemical or physical agents[END_REF][START_REF] Russell | Some aspects of thermal injury in Escherichia coli[END_REF][START_REF] Do | Damage to Escherichia coli on exposure to moist heat[END_REF][START_REF] Mc Allwood | Mechanism of thermal injury in Staphylococcus aureus: I. relationship between viability and leakage[END_REF][START_REF] Mc Allwood | Mechanisms of thermal injury in nonsporulating bacteria[END_REF][START_REF] Tomlins | Thermal injury and inactivation in vegetative bacteria. Inhibition and inactivation of vegetative microbes[END_REF][START_REF] Lr Beuchat | Injury and repair of gram-negative bacteria, with special consideration of the involvement of the cytoplasmic membrane[END_REF]. Such injury can be produced by physical processes such as extremely high temperatures (compared with the bacterial T OG ) and freezing [START_REF] Denver | Microbial susceptibility and resistance to chemical and physical agents[END_REF]. Membrane damage can be detected very readily by measuring the extent of intracellular material that leaks from the heated cells.

Overall, even if the rate and extent of leakage increases with the applied heating temperature, the correlation between loss of viability and membrane damage is poor [START_REF] Russell | Some aspects of thermal injury in Escherichia coli[END_REF][START_REF] Do | Damage to Escherichia coli on exposure to moist heat[END_REF][START_REF] Mc Allwood | Mechanism of thermal injury in Staphylococcus aureus: I. relationship between viability and leakage[END_REF][START_REF] Mc Allwood | Mechanisms of thermal injury in nonsporulating bacteria[END_REF]. For this reason, it seems that membrane damage is not the major cellular site responsible for inactivation. This is even more true for Gram-positive bacteria that are relatively more heat resistant than Gram-negative bacteria [START_REF] James | Taxonomy, role, and significance of microorganisms in foods[END_REF]. "Evaluation of heating effects on the morphology and membrane structure of Escherichia coli using electron paramagnetic resonance spectroscopy" (Biophysical Chemistry 2019) [START_REF] Tonyali | Evaluation of heating effects on the morphology and membrane structure of Escherichia coli using electron paramagnetic resonance spectroscopy[END_REF].

• Ribosome and rRNA degradation. Many studies have been performed on rRNA and ribosome to study their thermal stability in heated bacterial cells [START_REF] Mc Allwood | Mechanisms of thermal injury in nonsporulating bacteria[END_REF][START_REF] Tomlins | Thermal injury and inactivation in vegetative bacteria. Inhibition and inactivation of vegetative microbes[END_REF][START_REF] John | Repair of thermal injury of Staphylococcus aureus[END_REF][START_REF] Mc Allwood | Thermally induced ribonucleic acid degradation and leakage of substances from the metabolic pool in Staphylococcus aureus[END_REF][START_REF] Richard | Precursor ribosomal ribonucleic acid and ribosome accumulation in vivo during the recovery of Salmonella typhimurium from thermal injury[END_REF]. Mild heating produces degradation of rRNA [START_REF] Mc Allwood | Thermally induced ribonucleic acid degradation and leakage of substances from the metabolic pool in Staphylococcus aureus[END_REF], and the degradation of 30S ribosomal subunits seems to be very dependent on the concentration of salt ions in the cytoplasm [START_REF] Hurst | Stability of ribosomes of Staphylococcus aureus S6 sublethally heated in different buffers[END_REF][START_REF] Hurst | Revival of vegetative bacteria after sublethal heating[END_REF].

However, the rRNA degradation occurs before the loss of the cellular viability [START_REF] Ne Welker | Microbial endurance and resistance to heat stress[END_REF]. On the other hand, Tomlins & Ordal [START_REF] Tomlins | Thermal injury and inactivation in vegetative bacteria. Inhibition and inactivation of vegetative microbes[END_REF] found that the degradation of ribosomal RNA is reversible. Therefore, these are not considered to be primary causes of heat inactivation.

• DNA damage. Due to its thermal resistance, the denaturation of DNA can be considered to be a minor cause of the cell inactivation due to heat shocks [START_REF] Ne Welker | Microbial endurance and resistance to heat stress[END_REF].

However, there is a relationship between the bacterial sensitivities to ionizing radiation and mild heat shock [START_REF] Bridges | Correlation of bacterial sensitivities to ionizing radiation and mild heating[END_REF][START_REF] Pauling | Role of DNA ligase in the repair of single strand breaks induced in DNA by mild heating of Escherichia coli[END_REF][START_REF] John R Battista | The stress responses of Deinococcus radiodurans[END_REF].

• Protein denaturation. Structural and functional proteins are a prime target for heat inactivation. Protein denaturation occurs when cells are thermally stressed, with coagulation occurring at very high temperatures [START_REF] Mc Allwood | Mechanisms of thermal injury in nonsporulating bacteria[END_REF][START_REF] Harries | A note on some changes in the physical properties of Escherichia coli after heat treatment[END_REF]. There is a numerical correlation between the thermodynamic parameters of protein denaturation and the observed death rates of various organisms, suggesting a likely cause of cell death in mesophilic bacteria [START_REF] Rosenberg | Quantitative evidence for protein denaturation as the cause of thermal death[END_REF].

This relation is strengthened by the observation that, in thermophiles and hyperthermophiles, enzymes and proteins have an enhanced thermal stability and appear to function optimally at high temperatures. The heat stability of proteins from thermophiles is increased by the presence of a number of salt bridges and by the densely packed hydrophobic interior nature of the proteins 1 .

In conclusion, lethal effects of high temperatures are due to a combination of processes that involve simultaneously different parts of the bacteria. All the major cellular components are impacted, namely the outer cell barrier (especially for Gram-negative bacteria), ribosomal RNA, proteins and DNA. To react to thermal shock, the cells have different active mechanisms, in particular, there is a rapid induction of intracellular heat-shock proteins (HSPs). HSPs behave as molecular chaperones binding to the proteins that are becoming unstable due to the increase in temperature, helping to stabilize their structure and preventing their unfolding [START_REF] Georgopoulos | Properties of heat shock proteins of Escherichia coli and autoregulation of the heat shock response[END_REF][START_REF] Chen | Thermosensitivity of growth is determined by chaperonemediated proteome reallocation[END_REF] -see Figure 1 Moreover, extracellular alarmones may be produced to warn organisms of imminent inactivation. Exposure of cultures to increasing temperatures may enable the cells to adapt to higher, normally rapidly lethal temperatures [START_REF] Russel | Lethal Effects of Heat on Bacterial Physiology and Structure[END_REF].

The role of proteins: current situation

As described in the previous section, the proteome's thermal sensitivity has to play a key role as a determinant for most of the temperature-dependent whole-organism activities.

Different pictures have been proposed to link the degradation of the proteome to the upper limit of the cellular thermal niche, i.e. the cell's death temperature T CD .A first essential aspect is to quantify the proteome's thermal stability [START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF][START_REF] Dill | Physical limits of cells and proteomes[END_REF][START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF]. On one hand a proposed theoretical model [START_REF] Ghosh | Cellular proteomes have broad distributions of protein stability[END_REF][START_REF] Dill | Physical limits of cells and proteomes[END_REF] finds that cell death is linked to a global catastrophe of the proteome with proteins unfolding in a narrow range of temperatures near the T CD . This picture has being challenged recently by experimental investigations of E. coli lysates and cells, and based on different techniques such as limited proteolysis [START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF] or thermal proteome profile [START_REF] Mateus | Thermal Proteome Profiling in Bacteria: Probing Protein State in vivo[END_REF] combined with mass spectroscopy. According to these studies only a small set of proteins indeed unfolds at the cell death. Thermal adaptation would result from the preferential stabilization of a homologous subset of proteins, thus indicating that the heat sensitivity of cells can be explained by a small number of proteins that serve critical physiological roles.

Actually, the proteome's thermal stability is not the only physical determinant of the cell's growth rate, which is expected to depend on the rate of protein diffusion throughout the cell, the latter being often the limiting factor of the rates of cellular biochemical processes [START_REF] Paul E Schavemaker | How important is protein diffusion in prokaryotes[END_REF]. Protein diffusion depends in turn on the temperature, especially through the contribution of the intrinsic viscosity in the high-temperature range when biomolecules start to unfold. To date, the relationship between the diffusive dynamics of proteome and the thermal sensitivity of a cell has not yet been investigated, also due to the extremely difficult challenge to represent the motions of proteins in a crowded milieu cell's cytoplasm where local concentration may vary from 200 g/L up to 400 g/L [START_REF] John | Macromolecular crowding: an important but neglected aspect of the intracellular environment[END_REF]. Here, the protein diffusive dynamics is affected by several factors, such as the presence of steric barriers given by the other macromolecules, hydrodynamic and attractive interactions and spatial heterogeneity.

Our approach

To address this problem, we combined neutron scattering spectroscopy and multiscale molecular dynamics simulations to characterized the dynamical state of the E. coli's proteome in the thermal range between 276K and 360K, at increasing and decreasing temperatures to test its reversibility. We focus our studies on E. coli because they are the most studied bacteria, and many useful information are available in the literature that can be used to interpret the results. Moreover, this bacteria represent a good model for all bacterial cells, and, since we are mainly interested in the bulk properties of the system, i.e the dynamical behavior of the cytoplasmic proteins which are the most abundant in the cell, E. coli bacteria are a good reference also for the cytosol of eukaryotic cells.

Concerning the experimental techniques, in chapter 2, we presented the basics concepts on neutron scattering: introduction to neutron scattering theory, description of the model used for the interpretation of the data, and data reduction. Chapter 3 is dedicated to molecular dynamics simulation where we discussed about: the general aspects of molecular dynamics, the simulation of proteins, and coarse-grained simulations. Finally, in chapter 4 we report our new results and we discuss their interpretations and implications.

In the appendix, we reported also other two original results that we obtained, during the PhD, studying the effects of protein-ligand complexation on protein dynamics.

Chapter 2 Neutron Scattering

The aim of this chapter is to provide a short introduction to scattering theory. The main equations relevant for the scattering techniques used within this thesis are derived to provide a theoretical background. For more detailed explanations, literature e.g. by Lovesey [START_REF] Stephen W Lovesey | Theory of neutron scattering from condensed matter[END_REF], Bée [START_REF] Bée | Quasielastic neutron scattering[END_REF], Squires [START_REF] Squires | Introduction to the theory of thermal neutron scattering[END_REF], Boothroyd [START_REF] Andrew | Principles of Neutron Scattering from Condensed Matter[END_REF], can be consulted.

General concepts

The neutron is an electrically neutral particle with mass m n ⇡ 939.6MeV/c 2 , spin s = h/2 and a magnetic dipole moment µ n = 1.913µ N , where h = h 2π ⇡ 6.582µeV ns is the reduced Planck constant and µ N ⇡ 3.1525•10 2 µeV T 1 is the nuclear magneton. Neutrons are found in the atomic nuclei together with protons, where they are stable. In order to use them as a probe for neutron scattering experiments they have to be extracted from the nuclei. Free neutrons are not stable, they decay via a β decay process to a proton p + by emitting an electron e and an electron antineutrino ν e : n ! p + + e + ν e with a mean lifetime of approximately 880 seconds. Therefore, in order to perform experiments, they have to be produced steadily from a source that, depending on the energy of neutrons required for the measurements, should not be too far from the sample. To date, there are only two ways of producing high fluxes of free neutrons for experiments: through nuclear fission in nuclear reactors or by nuclear spallation in spallation sources.

The characteristics of a neutron flux produced by reactor based sources and spallation sources are different. A reactor produces neutrons at a constant rate and thus the flux of neutrons has no explicit time structure -these are called continuous or steady state neutron sources. Typical examples are the Institut Laue Langevin (ILL) reactor in France or the Forschungs-Neutronenquelle Heinz Maier-Leibnitz (FRM II) reactor in Germany. In contrast, spallation sources generally work with pulsed neutron beams and are thus called pulsed spallation sources. Typical examples are the ISIS in UK, the Spallation Neutron Source (SNS) in USA, and in the near future the European Spallation Source (ESS) in Sweden.

The experiments describes in this thesis were all carried out at the ILL, therefore we will briefly focus on nuclear fission being the production method employed at this facility. Nuclear fission, specifically neutron-induced fission in a nuclear reactor is the process by which a heavy nucleus like 235 U absorbs a neutron and subsequently splits into two or more lighter nuclei accompanied by the release, on average, of 2 or 3 neutrons per fission with typical kinetic energies of about 2 MeV per neutron (Fig. 2.1). These neutrons are released after the fission because the neutron/proton ratio for stable nuclei increases with the increase of their atomic number Z and, as consequence, when heavy nucleus splits into stable lighter nuclei there is an excess of neutrons that will be released with the fission. The core of a nuclear reactor contains the fissile fuel elements and it is settled so that it there is an high probability that neutrons released after the fission of a nucleus induce at least one additional fission by collision with other nuclei. This in turn induces another fission, causing a chain reaction.

In the case of nuclear reactors designed for experiments, like the one of ILL, the core is built so that the majority of the excess neutrons which do not participate in fission reactions are extracted from the reactor through special guides, usually coated with nickel and titanium, in order to guide the neutrons from the source to the instruments. Moreover, since the energies of the neutrons generally needed for scattering experiments are usually lower than few eV, to be suitable for experiments, the neutrons extracted from the reactor must be slowed down to reduce their kinetic energy. This is achieved with moderators, which are mediums that slow the neutrons down by repeated collisions with their nuclei. After many collisions the neutrons reach thermal equilibrium with the moderator and the flux per unit energy φ (E) of out-coming neutrons can described using a Maxwellian distribution, whose maximum is given by the moderator temperature T M [START_REF] Andrew | Principles of Neutron Scattering from Condensed Matter[END_REF]:

φ (E) ∝ E • exp ✓ E k B T M ◆ (2.1)
and the mean energy is determined by the temperature of the moderator. Examples of actual moderator substances are liquid hydrogen at 25K, liquid water at 300K, and solid graphite at 2,400K1 . For obvious reasons, neutrons emerging from these moderators are termed cold, thermal, and hot, respectively. In these energy ranges, neutrons can be treated as non-relativistic particles. Thus, given the wave-particle duality, the kinetic energy E of a neutron can either be described as a particle by its mass m n and momentum hk or as a wave with a wavelength λ .

E = h2 k 2 2m n = h 2 2m n λ 2 (2.2)
In particular, thermal neutrons (e.g. moderated by a moderator with T M = 293K) have a mean velocity of v ⇡ 2.20km/s and a corresponding wavelength of λ ⇡ 1.8 Å -see Tab. 2.1. Since this wavelength is comparable to the interatomic distances of liquids and solids, thermal neutrons provide an ideal tool for probing the microscopic properties of these materials. In addition, the energy of thermal neutrons is comparable to thermal excitations and can therefore be used to probe molecular vibrations, lattice excitations as well as atomic dynamics.

Combining the average neutron speed with the mean lifetime of the neutron, it is clear that the construction of instruments situated several hundred meters away from the neutron source is feasible even for cold neutrons. This offers the possibility to construct many instruments using one neutron source and to use e.g. curved neutron guides blocking the direct view from the instrument towards the reactor core and therefore reducing the background signal.

A diagram showing the key components of the research reactor at the ILL is provided in Fig. 2.3. The reactor operates at an output power of 56MW and produces a steady flux of neutrons in the moderator region of 1.5 • 10 15 cm 2 s 1 . There are a number of cold, thermal, and hot moderators, and the neutrons from these are delivered to over forty instruments designed for a wide range of scattering and fundamental physics experiments.

Neutron Scattering -Basics

In neutron scattering experiments a beam of neutrons with energy E i is used to interact with a sample to obtain information on the dynamical structure and the motions of its nuclei. To this end, the flux Φ of incident neutrons is measured as well as the rate J of the scattered neutrons with a final energy between E f and E f + ∆E f that are revealed, with a certain efficiency η, by a detector placed far from the point of interaction that covers a small solid angle ∆Ω in the direction (θ , φ ). In general, the scattering is axial symmetric i.e. the rate J does not depend on the φ direction, and it can be described as follow

J(θ , E i , E f )=η(θ ) • Φ(E i )  d 2 σ dΩ dE f (θ , E i , E f ) • ∆E f • ∆Ω (2.3)
where d2 σ dΩ dE f is the double differential cross-section, which represents the probability that a neutron in the initial state |ψ i i with an energy E i and wave-vector k i passes over into the final state ψ f ↵ characterized with an energy E f and the wave-vector k f that forms an angle θ with k i . The probability of the transition from |ψ i i to ψ f ↵ can be described by Fermi's Golden Rule.

d 2 σ dΩ dE f = m 2π h2 k f k i ⌦ ψ f V N (r) |ψ i i 2 δ (E + E i E f ) (2.4)
where δ (E + E i E f ) represents the energy conservation (E is the energy transferred from the neutron to the sample and it should be always equal to E i E f ) and

k f k i is nec- essary for the normalization of the neutron flux. V N is the so called pseudo potential with V N (r)= 2π h2 m ∑ α b α δ (r rα )
where N is the number of nuclei in the sample, rα is the position of the α th particle, and b α is the corresponding scattering length which represents the strength of the potential 2 . This pseudo potential causes the same scattering as the real potential, but is weak enough to be treated by perturbation expansion derived by Max Born. In first approximation, the Born expansion states that if |ψ i i can be described by a plane wave with wave vector k i , the final state ψ f ↵ is as well a plane wave with wave vector k f , and for the cross-section we obtain the following relation

d 2 σ dΩ dE f = k f k i * ∑ α, β b ⇤ α e i q•r α • b β e i q•r β + δ (E + E i E f ) (2.5)
where q = ki k f is the wave-vector transferred from the neutron to the sample. Now, take into account the integral form of the δ -function:

δ (E + E i E f )= 1 2π h Z +∞ ∞ dt • exp ✓ it • E + E i E f h ◆ (2.6)
it can be substituted in eq. (2.5) leading to [START_REF] Squires | Introduction to the theory of thermal neutron scattering[END_REF]:

d 2 σ dΩ dE f = 1 2π h k f k i Z +∞ ∞ dt • e itE/h ∑ α, β D b ⇤ α b β e i q•r α (t) • e i q•r β (0) E (2.7)
It is important to observe that q and rα are actually operators and that generally they do not commute. As a consequence, [r α , rα (t)] 6 = 0 and therefore the product of e i q•r α (0) and e i q•r α (t) is not equal to e i q•(r α (0) r α (t)) .

In the classical limit wave-vectors and the atomic positions are not operators, which means that we replace q ! q and ri ! r i , and consequently:

e iq•r α (t) e iq•r β (0) = e iq•[r α (t) r β (0)] (2.8)
Therefore for the cross-section we obtain

d 2 σ dΩ dE f = 1 2π h k f k i Z +∞ ∞ dt • e itE/h ∑ α, β D b ⇤ α b β e iq•[r α (t) r β (0)] E cl (2.9)
where h...i cl means that a classical average is taken. In the following, the index "cl" will be omitted.

Coherent and Incoherent Scattering

The spins of the nuclei in the sample and the spins of the neutrons in the beam are uncorrelated in our experiments. Assuming two scattering atoms α, β of the same isotopes having the scattering lengths b α and b β due to their spin states. The evaluation of ⌦ b ⇤ α b β ↵ leads to:

α 6 = β : ⌦ b ⇤ α b β ↵ = hb ⇤ α i ⌦ b β ↵ = hbi 2 α = β : ⌦ b ⇤ α b β ↵ = ⌦ b 2 ↵ (2.10)
Combining the two equations above,

⌦ b ⇤ α b β ↵ can be split in two parts, reading ⌦ b ⇤ α b β ↵ = hbi 2 + δ αβ (hbi 2 ⌦ b 2 ↵ )= 1 4π • (σ coh + δ αβ σ inc ) (2.11) with σ coh = 4π ⌦ b 2 ↵ and σ inc = 4π( ⌦ b 2 ↵ hbi 2
). With this result the scattering crosssection in eq. (2.9) can be split into two components:

d 2 σ dΩ dE f = N 4π k f k i [σ coh S coh (q, E)+σ inc S inc (q, E)] (2.12)
where N is the total number of nuclei in the scattering system and

S coh (q, E)= 1 2π h Z +∞ ∞ dt • e i tE h 1 N ∑ α, β D e iq•[r α (t) r β (0)] E (2.13) S inc (q, E)= 1 2π h Z +∞ ∞ dt e i tE h 1 N ∑ α D e iq•[r α (t) r α (0)] E (2.14)
S coh (q, E) and S inc (q, E) are respectively the coherent and incoherent dynamic structure factor (also known as scattering functions). Therefore, the coherent scattering depends on the self-correlation between the positions of the same nucleus at different times and on the correlation between the positions of different nuclei at different times, leading to interference effects. In contrast, the incoherent scattering depends only on the selfcorrelation between the positions of the same nucleus at different times which do not give interference effects. Thus, coherent scattering can give information on structure and collective motions (correlations between all scatterers), whereas incoherent scattering gives information about the space evolution in time (self-correlation), leading to a probe of the local dynamics of the sample. The scattering function S coh (q, E) and S inc (q, E) can be also expressed as the Fourier transform in time of the coherent and the incoherent intermediate scattering function I(q,t) defined as

I coh (q,t)= 1 N ∑ α, β D e iq•[r α (t) r β (0)] E (2.15) I inc (q,t)= 1 N ∑ α D e iq•[r α (t) r α (0)] E (2.16)
according to which:

S coh (q, E)= 1 2π h Z +∞ ∞ dt • I coh (q,t) (2.17) S inc (q, E)= 1 2π h Z +∞ ∞ dt • I inc (q,t) (2.18)
It is important to observe that, if σ inc ⌧ σ coh , then the incoherent scattering domi- nates the scattering signal and we can approximate

d 2 σ dΩ dE f ⇡ N 4π • k f k i • σ inc S inc (q, E) (2.19)
This is a good approximation for biological samples, like proteins or bacteria. As reported in Tab. 2.2, hydrogen H has a very large incoherent cross section of about 80 barns and in general, biological systems are approximately formed by half of H atoms. The remaining atoms mainly consist of C, N and O which have all very small incoherent cross sections (less than 0.5 barns), thus the majority of the signal comes from the incoherent scattering of H atoms in the sample. In most systems, the surrounding water (H2O) present in the sample will add an important contribution to the signal, thus often heavy water D2O is used as medium since its incoherent scattering cross section is 40 times smaller with 2 barns. Deuterium can also be used as substitute of hydrogen to mask specific regions in the protein.

Therefore, in what follows, we will only consider the incoherent scattering and the index "inc" will be omitted. Finally, in many neutron scattering experiments it is useful to split the dynamical structure factor in three components, depending on the amount of energy transferred with the scattering:

1. Elastic incoherent neutron scattering (EINS): neutrons do not exchange energy with the sample, i.e. E ⇡ 0.

2. Quasi-elastic neutron scattering (QENS): only a small amount of energy is exchanged between the neutron and the sample (typically less then few meV), which produce a broadening of the elastic peak.

3. Inelastic neutron scattering (INS) which appears as additional peaks centered at E 6 = 0 well separated from the elastic peak that correspond to specific modes or excitations within the sample.

A schematic is shown in Figure 2.5.

Elastic Incoherent Neutron Scattering (EINS)

In the case of elastic incoherent neutron scattering (EINS), E ⇡ 0, i.e. the transferred energy is comprised between 0 and ∆E, where ∆E is the accuracy on the measurement of the energy. As it will be described with more detail in section 2.3.2, the measured dynamical structure factor S exp depends on the theoretical scattering function S th (which comes from the sample-neutron interaction), and the resolution function of the instrument as follow:

S exp (q, E = 0; ∆ω)=S th (q, 0) ⌦ R(q, E; ∆E)= Z +∞ ∞ dE 0 S th (q, E 0 ) R(q, E 0 ; ∆E) (2.20)
and, considering that in our experiments the resolution of the instruments is generally well represented by a Gaussian where ∆E corresponds to its Full Width Half Maximum (FWHM), the measured scattering function S exp (q, E = 0; ∆E) is approximately equal to the intermediate scattering function I exp (q,t R ) [START_REF] Doster | TOF-elastic resolution spectroscopy: time domain analysis of weakly scattering (biological) samples[END_REF]:

S exp (q, E = 0; ∆E) ⇡ I exp (q,t R )= 1 N N ∑ α=1 he iq•[r α (t R ) r α (0)] i (2.21)
where t R = 4h p π 1 ln 2/∆E is the time window fixed by the instrumental resolution and determines the experimentally observable motions [START_REF] Zorn | On the evaluation of neutron scattering elastic scan data[END_REF]. Therefore, in elastic neutron scattering the rate of scattered neutrons is proportional to the self intermediate scattering function I(q,t) and, as a consequence, we can focus directly on this function instead of the S exp (q, E = 0; ∆E).

If we define the displacement ∆r α (t R )=r α (0) r α (t R ), the self intermediate scattering function can be describe as:

I(q,t R )= 1 N N ∑ α=1 he iq•∆r α (t R ) i (2.22)
The index "exp" will be omitted in the following sections.

Gaussian Approximation (GA)

Starting from eq. 2.22, we can use the cumulant expansion for each term of the sum [START_REF] Squires | Introduction to the theory of thermal neutron scattering[END_REF]: .23) with:

he iq•∆r α (t R ) i = e q 2 • ρ 2 (t R )+q 4 • ρ 4 (t R ) ± ... ( 2
ρ 2 (t)= 1 2! hn q • ∆r 2 α (t)i ρ 4 (t)= 1 4! ⇥ hn q • ∆r 4 α (t)i 3 hn q • ∆r 2 α (t)i 2 ⇤ (2.24)
where n q is the direction of q and hn q • ∆r n α (t)i is the mean value of the n-th power of displacements in the direction of q.

Then, neglecting the terms of order higher than q 2 and assuming that, on average, only isotropic displacements take place3 (such that hn q • ∆r n α i = h∆r 2 α i/3), the intermediate scattering function can be written as:

I(q,t R ) ⇡ 1 N N ∑ α=1 e 1 6 q 2 •h∆r 2 α (t R )i (2.25)
Now, if we neglect also the motional heterogeneity of the H atoms in the system and we assume that the h∆r 2 α i of all the H are the same, i.e. the harmonic potential is equal for all the H atoms, the h∆r 2 α i = h∆r 2 i and we obtain the well known Gaussian Approximation (GA):

I(q,t R ) ⇡ I 0 • e 1 6 q 2 •h∆r 2 (t R )i (2.26)
where I 0 is the value of I(q,t R ) for |q| = 0 and

ln(S) ⇡ 1 6 q 2 h∆r 2 (t R )i + ln(I 0 ) (2.27)
It is important to observe that the GA is generally only valid in a restricted region of q (specifically at low q values) since it neglects any effects from anharmonicity, heterogeneity or anisotropy. Each of these effects alone or any combination can lead to a divergence from the Gaussian behaviour. As shown by Tokuhisa et al. [START_REF] Tokuhisa | Non-Gaussian behavior of elastic incoherent neutron scattering profiles of proteins studied by molecular dynamics simulation[END_REF] and Vural et al. [START_REF] Vural | Motional displacements in proteins: The origin of wave-vector-dependent values[END_REF], usually the heterogeneity of motions is the biggest contribution to the non Gaussianity of the self intermediate scattering function [START_REF] Zeller | Analysis of elastic incoherent neutron scattering data beyond the Gaussian approximation[END_REF].

Non-Gaussian Scattering

This model tries to correct, separately, the non-Gaussian behavior of the self intermediate scattering function due to single-atom non-Gaussian scattering (anharmonic displacements) and dynamic heterogeneity [START_REF] Yi | Derivation of mean-square displacements for protein dynamics from elastic incoherent neutron scattering[END_REF].

Non-Gaussian single-atom dynamics leads to the contribution of higher-order terms in eq. (2.23) being non-negligible, that can be expanded as [START_REF] Yi | Derivation of mean-square displacements for protein dynamics from elastic incoherent neutron scattering[END_REF]:

he iq•∆r α (t R ) i = e 1 6 q 2 •h∆r 2 (t R )i " 1 + ∞ ∑ k=2 ρ k ( q 2 ) k # (2.28)
where we did the approximation of isotropic displacements. Now, retaining only the k = 2 term and neglect the motional heterogeneity, with eq. (2.24), I(q,t R ) can be written as [START_REF] Yi | Derivation of mean-square displacements for protein dynamics from elastic incoherent neutron scattering[END_REF]:

I(q,t R )=e 1 6 q 2 •h∆r 2 (t R )i  1 + h∆r 4 (t R )i h∆r 2 (t R )i 2 72 • q 4 (2.29)
On the other hand, the dynamical heterogeneity also leads to non-Gaussian behavior and, starting from eq. (2.25), the correction can be made as follows [START_REF] Yi | Derivation of mean-square displacements for protein dynamics from elastic incoherent neutron scattering[END_REF]:

I(q,t R ) ⇡ e 1 6 q 2 •h∆r 2 (t R )i • 1 N ∑ N α=1 e 1 6 q 2 • (h∆r 2 α (t R )i h∆r 2 (t R )i) = e 1 6 q 2 •h∆r 2 (t R )i • 1 N ∑ N α=1 ∑ ∞ α=1 µ(m) m! ⇣ q 2 6 ⌘ m ⇡ e 1 6 q 2 •h∆r 2 (t R )i • ⇣ 1 + q 4 72 σ 2 ⌘ (2.30)
where µ(m) is the m-th central moment of the distribution of h∆r 2 (t R )i and σ 2 is the variance:

σ 2 = 1 N ∑ α=1 h∆r 2 α (t R )i h∆r 2 (t R )i 2 (2.31)
and the last approximation in eq. (2.29) is valid if ( q 2 /6) m µ(m) ⌧ 1. Both eq. (2.29) and eq. (2.30) have similar expressions, indicating the general validity of the proposed q 4 model. Moreover, in systems in which dynamical heterogeneity is the dominant contribution to non-Gaussian behavior, the elastic scattering can in principle be used to obtain experimentally the variance and higher statistical moments of the distribution of the mean-square displacements of individual atoms. It should be noted that this q 4 correction to the GA does not presume any particular form of the distribution of MSD [START_REF] Yi | Derivation of mean-square displacements for protein dynamics from elastic incoherent neutron scattering[END_REF].

However, since it does not contain higher order terms further than q 4 , its applicability is limited in q [60].

Dynamical Heterogeneity -Distribution function approaches

In order to incorporate dynamical heterogeneity, the distribution function should be explicitly included for the calculation of the self intermediate scattering function [START_REF] Nakagawa | Protein dynamical heterogeneity derived from neutron incoherent elastic scattering[END_REF]:

I(q,t R )= Z ∞ 0 g(s) • e sq 2 ds (2.32)
where g(s) is the distribution function of MSD and s is h∆r 2 (t R )i/6. It is useful to underline that, with this approach, we have not a single value for the MSD, but a distribution of values. In particular, the mean value of the MSD, h∆r 2 (t R )i, is given by:

s = Z ∞ 0 s • g(s) ds (2.33)
where the h∆r 2 (t R )i = 6 • s.

Peters and Kneller Model: This model tries to correct the GA model describing the dynamical heterogeneity with a Gamma distribution for the individual MSD. Following Peters and Kneller [START_REF] Peters | Motional heterogeneity in human acetylcholinesterase revealed by a non-Gaussian model for elastic incoherent neutron scattering[END_REF], we define the dimensionless momentum transfer q = lq where l > 0 is a scale variable with the dimension of a length, and substitute s by λ = s/l 2 . Thus, from eq. 2.33, we obtain:

I(q,t R )= Z ∞ 0 g(λ ) • e λ q 2 ds = Z ∞ 0 β (βλ) β 1 Γ(β ) • e βλ • e λ q 2 dλ (2.34)
where Γ(β ) is the Gamma function and β is a parameter of the distribution such that 0 < β < 1. This integral can be solved analytically and yields a simple analytical form for the model I(q,t R ):

I(q,t R )= 1 ⇣ 1 + h∆r 2 (t R )i 2 q 2 6 β ⌘ β (2.

35)

Bimodal: This model tries to correct the GA model describing the dynamical heterogeneity by a bimodal distribution for the individual MSD [START_REF] Nakagawa | Protein dynamical heterogeneity derived from neutron incoherent elastic scattering[END_REF]:

g(s)=p • δ (s s F )+(1 p) • δ (s s R ) (2.36)
where s F = h∆r 2 F (t R )i/6 and s R = h∆r 2 R (t R )i/6 are two MSD and p is a ratio of the two contributions such that 0  p  1. The interpretation of this model is that the system can be roughly separated in two classes of structures, flexible (F) and rigid (R), such that h∆r 2 F (t R )i < h∆r 2 R (t R )i. This leads to:

I(q,t R )=I 0 • h p • e 1 6 q 2 h∆r 2 F i +(1 p) • e 1 6 q 2 h∆r 2 R i i (2.37)
or, simplifying the eq. (2.37) in order to make easier the fit, we have: [START_REF] Pauling | Role of DNA ligase in the repair of single strand breaks induced in DNA by mild heating of Escherichia coli[END_REF] where:

I(q,t R )=I F • e 1 6 q 2 h∆r 2 F i + I R • e 1 6 q 2 h∆r 2 R i (2.
p = I F I R + I F and I 0 = I R + I F (2.39)
The mean value of the MSD is then:

h∆r 2 i = p •h∆r 2 F i +(1 p) •h∆r 2 R i (2.40)
Two q-Ranges: It is clear that eq. (2.37) opens the door to another approximation. Indeed, if there exist two q-ranges, q F and q R , such that:

I F e 1 6 q 2 h∆r 2 F i >> I R e 1 6 q 2 h∆r 2 R i f or q 2{q} F I F e 1 6 q 2 h∆r 2 F i << I R e 1 6 q 2 h∆r 2 R i f or q 2{q} R (2.41)
the eq. (2.37) yields to:

S el ⇡ ( I F e 1 6 q 2 h∆r 2 F i f or q 2{q} F I R e 1 6 q 2 h∆r 2 R i f or q 2{q} R (2.42)
and in analogy with eq. (2.27), this results in two linear regimes:

ln(S el ) ⇡ ⇢ 1 6 q 2 h∆r 2 F i + ln(I F ) f or q 2{q} F 1 6 q 2 h∆r 2 R i + ln(I R ) f or q 2{q} R (2.43)
The interpretation of this model is that in the system there are two different dynamical processes which are visible at different time scales [START_REF] Lehnert | Thermal Motions in Bacteriorhodopsin at Different Hydration Levels Studied by Neutron Scattering: Correlation with Kinetics and Light-Induced Conformational Changes[END_REF][START_REF] Zanotti | Experimental evidence of a liquid-liquid transition in interfacial water[END_REF][START_REF] Zeller | Analysis of elastic incoherent neutron scattering data beyond the Gaussian approximation[END_REF].

Generalized Mean Squared Displacement

Hennig et al. in 2011 [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF] have further developed the Gaussian approximation by transferring the approach by Rahman et al. [START_REF] Azms Rahman | Theory of slow neutron scattering by liquids[END_REF] to a the elastic scattering function S(q, E = 0), where a general q-dependent mean square displacement

⌦ u 2 ↵ q is intro- duced by e 1 6 q 2 •hu 2 i q := S(q, E = 0) (2.44) thus ⌦ u 2 ↵ q = 6 • ln [S(q, E = 0)] q 2 (2.45)
We approximate ⌦ u 2 ↵ q by using a Taylor expansion around q = 0 up to the 3 rd order

⌦ u 2 ↵ q = 3 ∑ n=0 a n n! q n + O(q 4 ) (2.46)
with the coefficients

a n = lim q!0 d n dq n ⌦ u 2 ↵ q (2.47)
Then, using a model for the description of S(q, E ⇡ 0) (similar to what we will see in next section on QENS), it is possible to show that the coefficients a 1 and a 3 are zero:

⌦ u 2 ↵ q = a 0 + a 2 q 2 + O(q 4 ) (2.48)
therefore S(q, E = 0)=e 1 6 (a0 q 2 +a 2 q 4 +O(q 6 )) (2.49)

and, if we take the limit q ! 0, we obtain a mean square displacement ⌦ u 2 ↵ that does not depend on q:

⌦ u 2 ↵ = lim q!0 ⇢ 6 • ln [S(q, E = 0)] q 2 = a 0 (2.50)
From this result, it can be proved that it possible to decomposed ⌦ u 2 ↵ as follows [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF]: ↵ are two diffusive contributions: the first is due to the internal motions of the molecular subunit (e.g. conformational changes), the second arises from the roto-translation of the entire molecule described by apparent diffusion coefficient D app . In particular, ⌦ u 2 diff ↵ can be calculated from D app :

⌦ u 2 ↵ = ⌦ u 2 vib ↵ + ⌦ u 2 sub ↵ + ⌦ u 2 diff ↵ (2.
⌦ u 2 diff ↵ = 6t R D app (2.52)
and consequently, also the internal mean square displacement

⌦ u 2 int ↵ = ⌦ u 2 vib ↵ + ⌦ u 2 sub ↵ can be obtained as ⌦ u 2 int ↵ = ⌦ u 2 ↵ ⌦ u 2 diff ↵ (2.53)

Quasi-Elastic Neutron Scattering (QENS)

The section is based on references [START_REF] Bée | Quasielastic neutron scattering[END_REF][START_REF] Gerald R Kneller | Quasielastic Neutron Scattering[END_REF], if no other references are mentioned. The incoherent intermediate scattering function in eq. (2.16) can be written as

I(q,t)= 1 N ∑ α D e iq•R α (t) E (2.54)
where R α (t)=r α (t) r α (0) is the displacement of the α-th atom covered within the time t. For molecules, it is often possible to decompose the motion of an atom as a superposition of different diffusive modes acting on different time-scales. In particular, in the case of large molecules like proteins, we can rewrite R α (t) as a combination of translational ("trn"), rotational ("rot"), internal diffusive ("int") and vibrational ("vib") motions: R α (t)= r trn (t)+r rot (t)+r int (t)+r vib (t) (2.55) therefore, assuming that those contributions are independent, I(q,t) can then be decomposed likewise in a product with:

I(q,t)= I trn (q,t) • I rot (q,t) • I int (q,t) • I vib (q,t) (2.56)
and, consequently, for the incoherent scattering function holds relation:

S(q, E)= S trn (q, E) ⌦ S rot (q, E) ⌦ S int (q, E) ⌦ S vib (q, E) (2.57)
where ⌦ is the convolution with respect to E. This equation is the initial point for the evaluation of the QENS spectrum S(q, E). Each component in eq. 2.57 leads to a broadening of the quasielastic peak, which is illustrated by the orange curve in Fig. 2.5. As a consequence, the broadening of this peak gives information about the different types of diffusion within the probed system.

Translational Diffusion: In the simplest case, the component S trn (q, E) reflects the translational diffusion of the center of mass (CoM) of the molecule which is moving freely and can described by the Brownian motion. The probability of finding the CoM of the molecule at a position r trn (t) after a time t is given by the incoherent selfcorrelation function G(t,t):

G(r,t)= ( 4 π D t t) 3 2 • exp  r 2 4 D t t (2.58)
where D t is the translational diffusion coefficient. The intermediate scattering function is obtained by the Fourier transform in space leading to:

I(q,t)=e q 2 D t t (2.59)
The scattering function is obtained by the Fourier transform of eq. 2.59 in time which leads to a Lorentzian reading [START_REF] Grimaldo | Dynamics of proteins in solution[END_REF]:

S trn (q, E)= 1 π • γ t (q) E 2 + γ 2 t (q) (2.60)
where the parameter γ t (q) is the width of the Lorentzian and, in the case of free diffu- sion, it is equal to D t q 2 . It is possible to prove that, even for more complex types of diffusive processes (like jump-diffusion), the scattering function can be described with a single Lorentzian, but the width γ t (q) will have a different dependence on q (in the case of jump diffusion, γ t (q)= D t q 2 1+D t q 2 τ ) [START_REF] Grimaldo | Dynamics of proteins in solution[END_REF].

Global Dynamics: Global diffusion describes the superposition of translational and rotational diffusion . In 1999, Perez et al. [START_REF] Pérez | Evolution of the internal dynamics of two globular proteins from dry powder to solution[END_REF] showed that the translational and the rotational motions of proteins acts often on similar time-scales, therefore, the contribution to the dynamical structure factor for the rotational and translational diffusion of proteins cannot be separated. On the contrary, they can be approximated by a single Lorentzian:

S glb (q, E)= 1 π • γ g (q) E 2 + γ 2 g (q)
(2.61)

where the width γ g (q) depends on an apparent diffusion coefficient D app which, in turn, is a function of both the translational and the rotational diffusion coefficients [START_REF] Roosen-Runge | Protein selfdiffusion in crowded solutions[END_REF].

Internal Dynamics: Proteins in aqueous solution are not rigid since the peptide chain is able to move within a confined space. The corresponding structure factor can be approximated by

S int (q, E)=A 0 (q) • δ (E)+[1 A 0 (q)] L α (E) (2.62)
with the Kohlrausch-Williams-Watts function [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF]:

L α (E)= 1 2π Z +∞ ∞ dt • e itE/h • e |t Γ| α (2.63)
This function is a generalized Lorentzian function with a broad range of relaxation times occurring in the dynamics of the protein. 0  α  1 is a phenomenological parameter, and for α = 1, L 1 is a normal Lorentzian with width Γ. A 0 (q) is the so called elastic incoherent structure factor (EISF) which describes the confined internal dynamics [START_REF] Bée | Quasielastic neutron scattering[END_REF].

The EISF describes the internal dynamics of a protein more precisely. These internal dynamics are very complex and therefore difficult to evaluate analytically. Notwithstanding, it is possible to apply a simplified model in which the EISF A 0 (q) is assumed to be a simple sum of weighted distributions describing different motions that are not uncorrelated [START_REF] Dj Bicout | Incoherent neutron scattering functions for combined dynamics[END_REF]:

A 0 (q)=p +(1 p) ⇥ sA sph (q)+(1 s) A 3JD (q) ⇤ (2.64)
where p is the fraction of atoms which appears to be fixed within the time scale of the instrument, whereas s is the fraction of hydrogen atoms undergoing a three-site jump diffusion. A 3JD(q) is the intensity caused by the jump diffusion process of the H-atoms in methyl groups (-CH 3 ). Those atoms are located on a circle at an average distance a M ⇡ 1.715 Å from each other and they perform jumps of 120 around the 3-fold axis [START_REF] Bée | Quasielastic neutron scattering[END_REF][START_REF] Pérez | Evolution of the internal dynamics of two globular proteins from dry powder to solution[END_REF].

A 3JD (q, T )= 1 + 2 j 0 (q • a M ) 3 where j 0 (x)= sin x x (2.65)
The internal motions are additionally described by the model of an atom diffusing in an impermeable spherical volume of radius r. The corresponding scattering intensity A sph (q) is given by

A sph (q)= ✓ 3 • j 1 (q • r) q • r ◆ 2 where j 1 (x)= sin x x 2 cos x x (2.66)
Atomic Vibration: Atoms in the protein oscillate around their equilibrium position. The structure factor of these vibrations can be decomposed into an elastic and an in-elastic part

S vib (q, E)=e 1 6 hu 2 vib i•q 2 (2.67)
The inelastic component S inel (q, E) has pronounced peaks at distinct energies E, as illustrated by the green line in Fig. 2.5.

Data reduction

When performing QENS measurements, very large arrays containing the counts from every detector are recorded together with the time when each neutron was detected -the time information is essential to determine the energy of the neutrons. In order to obtain the scattering function described in the previous sections, before any data analysis can be performed, the data must be reduced to a treatable form: for every count collected from different detectors and time channels is associated with a transferred wave-vector q, and a transferred energy E. Then, the data are grouped and integrated with the desired binning of both energy and wave-vectors. More specifically, the Mantid routines [START_REF] Arnold | Mantid-Data analysis and visualization package for neutron scattering and µ SR experiments[END_REF] provided by the ILL facility were used to reduce data. The subsequent analysis and fitting was preformed by self-written python available on github: https://github.com/DanieleDiBari/NSAnalysis. In general, QENS spectra present two main issues: 1) there is a contribution to the measured signal coming from the aluminum sample holder; 2) instrument-dependent calibration accounting for different efficiency of different detectors and the instrument resolution.

Empty Cell Subtraction

Even if the scattering cross-sections of aluminum are quite small (see Table 2.2), the contribution to the raw spectra arises from the aluminum cell containing the sample is not negligible. In order to correct such a contribution from the total scattering function, the spectrum of the empty sample holder is measured in addition to the samples. However, when the gap between the cell walls is filled with the sample, the spectrum of the cell is slightly different from that of the empty sample holder. In fact, both the incident and the scattered beams are subject to absorption not only within the walls of the sample cell, but also within the sample. To account for this effect, if the atomic composition of the sample is known, it is possible to employ the Paalman-Pings correction [START_REF] Hh Paalman | Numerical evaluation of X-ray absorption factors for cylindrical samples and annular sample cells[END_REF]. In the following, the superscript refers to the scattering object, while the subscript indicates the absorbing object. In both the super-and the subscript, "c" indicates the cell, and "s" the sample. Therefore, the scattering intensity of the bare sample I s , without the contribution of the empty cell is given by

I s (q, E)=α sc • I sc sc (q, E) β sc • I c c (q, E) (2.68)
where I sc sc indicates the intensity after scattering and absorption from both the sample and the cell, meanwhile I c c is the scattering intensity of the empty cell. Concerning the parameters α sc and β sc , they are defined as follows [START_REF] Hh Paalman | Numerical evaluation of X-ray absorption factors for cylindrical samples and annular sample cells[END_REF]:

α sc = 1 A s sc β sc = 1 A s sc • A c sc A c c (2.69)
where A s sc , A c sc , and A c c are the absorption factors, also known as Paalman-Pings coefficients, which depend on the scattering angle θ , the transferred energy E, the atomic composition of the sample, and the shape of cell.

Nevertheless, for complex samples, like bacteria or other cells, the atomic composition could be difficult to estimate. In these cases, measuring the transmission T sc of the sample and the cell, it is possible to approximate eq. (2.68) as follow:

I s (q,t) ⇡ I sc sc (q,t) T sc • I c c (q,t) (2.70)

Calibration: detector efficiency and energy resolution

The measured intensity of the QENS data depends, practically, on the efficiency of each detector as shown in eq. (2.3). In general we can say that the recorded intensity S exp (q, E) can be described by: S exp (q, E)=η(q) • S th (q, E) (2.71) where S th (q, E) is the scattering function theoretically produced by the interaction between the neutron beam and the sample. In the ideal case of perfect detector efficiency η(q), S exp (q, E)=S th (q, E), therefore η(q)=1. For incoherent QENS, η(q) can be estimated by integrating S V exp (q, E) of Vanadium (V):

η(q)= Z +E 0 E 0 S V exp (q, E) • dE (2.72)
whit ( E 0 , E 0 ) that is the energy range used for the integration elastic peak of the Vanadium. Vanadium is chosen because of its large incoherent neutron scattering crosssection, and its purely elastic peak (modulated by the Debye-Waller factor, see below) over the observable energy transfer range, meaning that ideally all the incoherent scattering intensity is detected. Therefore, all the measured spectra are corrected for the detector efficiency by the q-wise normalization:

S th (q, E)= S exp (q, E) R +E 0 E 0 S V exp (q, E) • dE (2.73)
However another crucial instrument-dependent feature that often affects every measurements, and that should be included in eq. (2.71), is the resolution function of the instrument. This function filters the signal produced by the scattering accordingly to the limitations of the data acquisition system of the instrument, and it depends generally on the transferred wave-vector and the time, R(q,t). In particular, it holds the following:

I exp (q,t)=η(q) • R(q,t) • I th (q,t) (2.74)
where I exp and I th are, respectively, the measured intermediate scattering function and the theoretical one. To compare this equation with eq. (2.71), we can take the Fourier transform (in time) of eq. (2.74) and we obtain

S exp (q, E)=η(q) • R(q, E) ⌦ S th (q, E) (2.75)
where R(q, E) Fourier transform of R(q,t).

As for the determination of the detector efficiency, measuring the scattering of Vanadium is ideal for determining the resolution function, since it is a strong elastic incoherent scatterer (Table 2.2) and serves as a standard to determine both the resolution function and the efficiency of the detectors of a neutron scattering instrument. The theoretical incoherent scattering function of vanadium in the µeV energy regime is [START_REF] Bée | Quasielastic neutron scattering[END_REF]:

S V th (q, E)=e 1 6 hu 2 i T •q 2 • δ (E) (2.76)
where ⌦ u 2 ↵ T is the temperature dependent mean-squared displacement, which for T = 296 K is ⌦ u 2 ↵ T =(6.7 ± 0.6) • 10 3 Å2 [START_REF] Kamal | Neutron incoherent elastic scattering study of the temperature dependence of the Debye-Waller exponent in vanadium[END_REF]. Hence, for the q-range of typical neutron scattering spectrometers, 0.2 Å 1  q  2 Å 1 , we can assume that the ⌦ u 2 ↵ T exponent in eq. (2.76) is negligible at ambient temperatures. Consequently, the vanadium signal can be approximated by a delta function with a q-independent peak intensity which is exactly what we need to characterize the detectors efficiency and their resolution function. Combining eq. (2.75) and (2.76) we obtain indeed the following relation:

S V exp (q, E)=η(q) • R(q, E) ⌦ S V th (q, E) = η(q) • R(q, E) ⌦ δ (E) = η(q) • R(q, E) (2.77)
Chapter 3

Molecular dynamics simulations

One technique employed in this thesis is Molecular Dynamics (MD) simulation. It is one of the most extended computer simulation techniques in material science and biophysics. In short, MD is used to calculate the time evolution of a classical many-body system by numerically integrating Newton's equations of motion. It is currently applied to a large spectrum of systems, from nanomaterials to biomolecules. This method provides two key features: it offers high spatial resolution and allows for the calculation of transport properties and relaxation constants since it does not disrupt the kinetics of the system. However, when applied to large systems, the sampling is limited by the current computational resources, e.g. the typical limit for a system of 100,000 atoms is the microsecond timescale [START_REF] Kalimeri | Are thermophilic proteins rigid or flexible? An in silico investigation[END_REF].

In the first section of this chapter, an introduction to the basic concepts of the MD simulation is presented. In the second section an insight on the atomistic MD simulation of proteins is shown -since proteins are complex systems they require specific techniques developed ad hoc. Finally, in the last section, the coarse graining strategy to simulate complex systems is presented.

Introduction to molecular dynamics simulations

MD simulations are one of the principal tools in the theoretical study of biomolecules providing information on the relative positions of molecules and atoms as a function of time and thus their dynamics. In fact, MD is a computational method which allows to calculate the time dependent behavior of complex molecular systems (e.g. proteins, nucleic acids, etc.) and hence it is generally used to describe such systems in terms of a realistic atomic model, with the aim to understand and predict macroscopic properties based on detailed knowledge on an atomic scale. Indeed, starting from an atomistic level, MD simulations are used to predict and better understand the properties of complex materials. In this way MD provide a direct route from the microscopic details of a system (the masses of the atoms, the interactions between them, etc.) to macroscopic properties of experimental interest (the equation of state, transport coefficients and so on). In particular, usually biomolecular MD simulations are used to gain insight into ligand binding, enzymatic activities, signalling mechanisms and protein folding [START_REF] Kukol | NAMD-VMD tutorial[END_REF]. Additionally simulations are valuable tools for the refinement of electron microscopy, x-ray, neutron scattering or other spectroscopic data in order to obtain more accurate molecular structures and used to interpret experimental results.

Actually, MD simulations compute the motions of individual molecules for a classical many-body system in order to describe the equilibrium and transport properties of solids, liquids and gasses. Although this modelling of the matter at the microscopic level must be, in principle, based on quantum mechanics, MD generally adopts a classical point of view. In this context, the word classical means that the nuclear motions of the constituent particles obey the laws of classical mechanics (the motions are described by the second Newton's law). This is an excellent approximation for a wide range of materials.

Hence, in MD neither relativistic nor quantum effects are considered:

Special relativity does not allow information to travel faster than light; MD simulations assume forces with an infinite speed of propagation.

Quantum mechanics has at its base the uncertainty principle; MD requires, and provides, complete information about position and momentum at all times.

In practice, the phenomena studied by MD simulations are those where relativistic effects are not observed and quantum effects can, if necessary, be incorporated as semiclassical corrections derived from quantum theory. 1MD simulations allow to calculate several properties of many-particle systems. However, not all properties can be directly measured in a simulation. Conversely, most of the quantities that can be measured in a simulation do not correspond to properties that are measured in real experiments. Actually, molecular simulations generate information at the microscopic level (atomic and molecular positions, velocities, etc.) and the conversion of this very detailed information into macroscopic terms (pressure, internal energy, etc.) is the field of the statistical mechanics. Thus, the language of statistical mechanics is necessary to use these simulations as the numerical counterpart of experiments.

In this context, it is useful to see that there is a direct connection between MD simulations and the microcanonical ensemble of statistical mechanics. Indeed, the microcanonical ensemble consists of all microscopic states (r N (t), p N (t)) on the constant energy hypersurface H(r N (t), p N (t)) = E. On the other hand in the classical Hamiltonian mechanics the equations of motion conserve the total energy:

dH dt = 0 =) H(r N (t), p N (t)) = const. (3.1)
This suggests a link between the microcanonical ensemble and classical Hamiltonian mechanics. For a system that evolves according to Hamilton's equations of motion, a trajectory computed with these equations, namely:

ṙi (t)= ∂ H ∂ p i ṗi (t)= ∂ H ∂ r i (3.2)
will generate microscopic configurations belonging to the microcanonical ensemble with the constant energy E. Thus, if after a long time a system with energy E is able to visit practically all the configurations on the constant energy hypersurface, the dynamical evolution of this system can be used to generate a microcanonical ensemble.

A system that has this property is said to be ergodic.

In general this dynamical approach, that is at the basis of MD simulations, provides a powerful method for generating an ensemble and its averages. 2 Thus MD simulations have evolved into one of the most widely used techniques for solving statistical mechanical problems. ). In an ergodic system, all the trajectories in this figure represent different sections of a single long trajectory. Indeed, if the system is ergodic, the single long trajectory would eventually pass through (or arbitrarily near) all states. A substantial region of cyclical trajectories, and a barrier region leading to bottleneck, are shaded. Source: Allen and Tildesley, Computer Simulation of Liquids (1st edition, 1987) [START_REF] Allen | Computer Simulation of Liquids[END_REF].

Given an ergodic trajectory, microcanonical phase space averages can be replaced by time averages over the trajectory according to:

hAi ⌘ R dx A(x) δ (H(x) E) R dx δ (H(x) E) = lim τ!∞ 1 τ Z τ 0 dt A[x(t)] ⌘ Ā (3.3)
where x(t) is a representative point of the phase space defined as x(t)=( r N (t), p N (t)). This formula can be discretized for MD simulations as follows:

hAi = 1 M M ∑ n=0 A(x n ∆t ) (3.4)
The discretization derives from the fact that the equations of motion are solved numerically using some numerical integrators that generate phase space vectors at discrete times that are multiples of a fundamental time discretization parameter ∆t, known as time step. Starting from x 0 , the vectors x n ∆t (where n = 0 ,...,M with M as the total number of integration steps) are generated by applying the integrator iteratively3 [START_REF] Di | Molecular Dynamics Simulation -A brief introduction of selected concepts[END_REF].

x(t) ! t ! n ∆t x n ∆t (3.5)
By the way, adopting the classical point of view, doing a MD simulation basically means solving Newton's equations of motion for a system of N interacting atoms:

m i ri = F i [r N (t)] i = 1 ,...,N (3.6)
where the forces are derived as the negative derivatives of the potential energy function U[r N (t)]:

F i [r N (t)] = ∇ r i U[r N (t)] (3.7)
The equations (3.6) and (3.7) are solved simultaneously at every time step. The system is supervised for a certain lapse of time, taking care that temperature and pressure remain at the required values and the coordinates are written to an output file at regular intervals. The coordinates, as a function of time, represent a trajectory of the system. After the initial changes, the system will usually reach an equilibrium state. By averaging over an equilibrium trajectory, many macroscopic properties can be extracted from an output file.

Molecular dynamics simulations of proteins

In this section several useful concepts are presented to outline how MD simulations of macromolecular systems generally work. Firstly, an introduction to the problem of the calculation of the forces is set: classical MD force field and boundary conditions, with a consideration of the special case of the long range Coulomb interactions, are shown. This is followed by a brief introduction to the numerical integration of (3.6) with mention to a simulation strategy for controlling temperature and pressure. Finally, the problem concerning the creation of the initial state is discussed: starting structures, solvation, minimization and equilibration.

Calculation of the forces

The potential energy is one of the most crucial parts of the simulation because it must faithfully represent the interaction between atoms, cast in the form of a simple mathematical function that can be calculated quickly. Indeed, the computation of the forces acting on every particle is the most time-consuming task of almost all MD simulations.

As biological systems involve many atoms of different types, a quantum mechanical treatment of these atoms is not feasible. The usual way to solve them is to use empirical potential energy functions, conventionally called force fields, which are computationally less expensive, but involve numerous approximations leading to certain limitations. 4 These functions and parameters have been derived from experimental results and quantum mechanical calculations of small model compounds. They are often refined by the use of computer simulations to compare calculated condensed phase properties with experiment. Current generation force fields provide a reasonable good compromise between accuracy and computational efficiency. Among the most commonly used potential energy functions are the AMBER, CHARMM, GROMOS and OPLS/AMBER force fields. One of the most important limitations of the empirical force fields is that no drastic changes in the electronic structure are allowed. i.e. no events like bond making or breaking can be modeled [START_REF] Gkeka | Molecular Dynamics simulations of lysozyme in water[END_REF].

The force field

As shown in eq. (3.7), the force acting on an atom i is calculated as the negative gradient of a scalar potential energy function U that depends on all atomic positions and, thereby, couples the motion of atoms. For systems of biomolecules, this potential energy function is usually divided into two parts:

U = U bonded + U non-bonded (3.8)
The bonded potential U bonded involves 2 , 3, and 4-body interactions of covalently bonded atoms, with O(N) terms in the summation. 5 The non-bonded potential U non-bonded involves long-range interactions between all pairs of atoms (usually excluding pairs of atoms already involved in a bonded term), with O(N 2 ) terms in the summation, although fast evaluation techniques are used to compute good approximations to their contributions to the potential with O(N) or O(N log N) computational cost. The different terms will be explained in more detail in the following sections.

Bonded potential terms

The bonded potential describes the stretching, bending, and torsional of the covalent bonds.

I Bond stretching:

The bond stretching term is a 2-body potential, generally assumed to be harmonic, that describes the vibrational motion between a pair of covalently bonded atoms:

U b = k b (b b 0 ) 2 (3.9)
where b is the distance between the two atoms. Two parameters characterize each bonded interaction: b 0 the average distance between them and a force constant k b .

I Angle bending:

The angle bending term describes the force originating from the deformation of the valence angles between three covalently bonded atoms (3-body interactions).

The angle bending term is described using a harmonic potential:

U θ = k θ (θ θ 0 ) 2 (3.10)
where θ is the angle between three atoms. Two parameters characterize each angle in the system: the reference angle θ 0 and a force constant k θ .

I Torsional terms:

The torsional terms are weaker than the bond stretching and angle bending terms. They describe the barriers to rotations existing between four bonded atoms (4body interaction). There are two types of torsional terms: proper and improper dihedrals. Proper torsional potentials are described by a cosine function:

U φ = k φ [1 + cos(n φ δ )] (3.11)
where φ is the angle between the planes formed by the first and the last three of the four atoms. Three parameters characterize this interaction: δ sets the minimum energy angle, k φ is a force constant and n is the periodicity.

The improper dihedral term is designed both to maintain chirality about a tetrahedral heavy atom and to maintain planarity about certain atoms. The potential is described by a harmonic function:

U ω = k ω (ω ω 0 ) 2 (3.12)
where ω is the angle between the plane formed by the central atom and two peripheral atoms and the plane formed by the peripheral atoms (see Fig. 3.2).

Non-bonded potential terms

The non-bonded potential describes the van der Waals forces and the electrostatic interactions between the atoms.

I Van der Waals interactions:

The van der Waals force acts on atoms in close proximity. It is strongly repulsive at short range and weakly attractive at medium range. The interaction is described by a Lennard-Jones potential:

U VdW = 4 ε  ⇣ σ r ⌘ 12 ⇣ σ r ⌘ 6 (3.13)
where r is the distance between two atoms. It is parametrized by σ : the collision parameter (the separation for which the energy is zero) and ε the depth of the potential well. The Lennard-Jones potential approaches 0 rapidly as r increases, so it is usually truncated (smoothly shifted) to 0 past a cutoff radius, requiring O(N) computational cost.

I Electrostatic interactions:

Finally, the long distance electrostatic interaction between two atoms is described by Coulomb's law:

U el = ε 1 4 • q 1 q 2 4 πε 0 r 12 (3.14)
where q 1 and q 2 are the charges of both atoms and r 12 the distance between them, while ε 0 is the electric susceptibility of vacuum. The parameter ε 1 4 is a unitless scaling factor whose value is 1, except for a modified 1 4 interaction, where the pair of atoms is separated by a sequence of three covalent bonds (so that the atoms might also be involved in a torsion angle interaction), in which case ε 1 4 = ε, for a fixed constant 0  ε  1. Although the electrostatic potential may be computed with a cutoff like the Lennard-Jones potential, the r 1 potential approaches 0 much more slowly than the r 6 potential, so neglecting the long range electrostatic terms can degrade qualitative results, especially for highly charged systems. There are other fast evaluation methods that approximate the contribution to the long range electrostatic terms that require O(N) or O(N log N) computational cost, depending on the method.

Potential energy function

Finally, the equation for the potential energy describing the force field can be expressed as:

U = ∑ bonds b k b (b b 0 ) 2 + ∑ angles θ k θ (θ θ 0 ) 2 + + ∑ proper dihedrals k φ [1 + cos(n φ δ )] + ∑ improper dihedrals k ω (ω ω 0 ) 2 + 9 > > > > = > > > > ; BONDED INTERACTIONS (3.15a) + ∑ i, j j<i 4 ε  ⇣ σ r ⌘ 12 ⇣ σ r ⌘ 6 + ∑ i, j j<i ε 1 4 • q 1 q 2 4 πε 0 r 12 =) NON-BONDED INTERACTIONS (3.15b)

Boundary condition

To avoid surface effects at the boundary of the simulated system, periodic boundary conditions are often used in MD simulations; the particles are enclosed in a cell that is replicated to infinity by periodic translations. A particle that leaves the cell on one side is replaced by a copy entering the cell on the opposite side, and each particle is subject to the potential from all other particles in the system including images in the surrounding cells, thus entirely eliminating surface effects (but not finite-size effects). Because every cell is an identical copy of all the others, all the image particles move together, consequently they should be represented only once inside the molecular dynamics code. However, because van der Waals and electrostatic interactions exist between every non-bonded pair of atoms in the system (including those in neighboring cells) computing the long-range interaction exactly is unfeasible. To perform this computation, the van der Waals interaction is spatially truncated at a user-specified cutoff distance. For a simulation using periodic boundary conditions, the system periodicity is exploited to compute the full (non-truncated) electrostatic interaction with minimal additional cost using the Particle-Mesh Ewald (PME) method described in the next paragraph.

Full Electrostatic Computation

Ewald summation is a description of the long-range electrostatic interactions for a spatially limited system with periodic boundary conditions. The infinite sum of chargecharge interactions for a charge-neutral system is conditionally convergent, meaning that the result of the summation depends on the order in which it is taken. Ewald summation specifies the order as follows: sum over each box first, then sum over spheres of boxes of increasingly larger radii. Ewald summation is considered more reliable than a cutoff scheme, although it is noted that the artificial periodicity can lead to bias in free energy, and can artificially stabilize a protein that should have unfolded quickly. The PME method is a fast numerical method to compute the Ewald sum. The cost of PME is proportional to N log N and the time reduction is significant even for a small system of several hundred atoms. The strict conservation of energy resulting from the computed force is crucial and is strongly assisted by maintaining the symplecticness of the integrator, as discussed further below. However the PME method does not conserve energy and momentum simultaneously, but momentum conservation can be enforced by subtracting the net force from the reciprocal sum computation, albeit at the cost of a small long-time energy drift.

Numerical Integration

Biomolecular simulations often require millions of time steps. Furthermore, biological systems are chaotic; trajectories starting from slightly different initial conditions diverge exponentially fast and after a few picoseconds are completely uncorrelated. However, highly accurate trajectories are not normally a goal for biomolecular simulations; more important is a proper sampling of phase space. Therefore, for constant energy (NVE ensemble) simulations, the key features of an integrator are not only how accurate it is locally, but also how efficient it is, and how well it preserves the fundamental dynamical properties, such as energy, momentum, time-reversibility, and symplecticness. The time evolution of a strict Hamiltonian system is symplectic. A consequence of this is the conservation of phase space volume along the trajectory, that is, the enforcement of the Liouville theorem. To a large extent, the trajectories computed by numerical integrators observing symplecticness represent the solution of a closely related problem that is still Hamiltonian. Because of this, the errors, unavoidably generated by an integrator at each time step, accumulate imperceptibly slowly, resulting in a very small long-time energy drift, if there is any at all. Artificial measures to conserve energy, for example, scaling the velocity at each time step so that the total energy is constant, lead to biased phase space sampling of the constant energy surface; in contrast, there has been no evidence that symplectic integrators have this problem.

One of the most used integratos for NVE simulations is the velocity Verlet [START_REF] Verlet | Computer" experiments" on classical fluids. i. thermodynamical properties of lennard-jones molecules[END_REF]. This method obtains the position and velocity at the next time step (r n+1 , v n+1 ) from the current one (r n , v n ), assuming the force F n = F(r n ) is already computed, in the following way:

half-kick ! v n+ 1 2 = v n + 0.5 • ∆t • F n /m drift ! r n+1 = r n + ∆t • v n+ 1 2 compute force ! F n+1 = F(r n+1 ) half-kick ! v n+1 = v n+ 1 2 + 0.5 • ∆t • F n+1 /m
where m is the mass. The Verlet method is symplectic and time reversible, conserves linear and angular momentum, and requires only one force evaluation for each time step. For a fixed time period, the method exhibits a (global) error proportional to ∆t 2 . More accurate (higher order) methods are desirable if they can increase the time step per force evaluation. However, higher order Runge-Kutta type methods, symplectic or not, are not suitable for biomolecular simulations because they require several force evaluations for each time step and force evaluation is by far the most time-consuming task in molecular dynamics simulations. Gear type predictor-corrector methods, or linear multistep methods in general, are not symplectic. Hence, no symplectic method has been found as yet that is both more accurate than the Verlet method and as practical for biomolecular simulations.

On the other hand, it is also possible to employ multiple-time-stepping methods to improve integration efficiency. Because the biomolecular interactions described in eq. ( 3.15) generally act on different time scales, this allows us to compute the slowervarying forces less frequently than faster ones. This idea can be implemented by three levels of integration loops. The inner loop uses only bonded forces to advance the system, the middle loop uses Lennard-Jones and short-range electrostatic forces, and the outer loop uses long-range electrostatic forces. 6Using this multiple time-step approach can increase computational efficiency by a factor of 2, however the longest time step is limited by resonance. 7

NVT Ensemble Simulations

A fundamental requirement for an integrator is to generate the correct ensemble distribution for the specified temperature and pressure in an appropriate way. For this purpose the Newtonian equations of motion (3.6) should be modified "mildly" so that the computed short-time trajectory can still be interpreted in a conventional way. To generate the correct ensemble distribution, the system is coupled to a reservoir, with the coupling being either deterministic or stochastic. Deterministic couplings generally have some conserved quantities (similar to total energy), the monitoring of which can provide some confidence in the simulation. NAMD uses a stochastic coupling approach because it is easier to implement and the friction terms tend to enhance the dynamical stability.

The (stochastic) Langevin equation is used in NAMD to generate the Boltzmann distribution for canonical (NVT) ensemble simulations. The generic Langevin equation is:

m i ri (t)=F i [r N (t)] γ ṙi (t)+ r 2 γ k B T m i R(t) (3.16)
where γ is the friction coefficient, k B is the Boltzmann constant, T is the temperature and R(t) is a univariate Gaussian random process. Coupling to the reservoir is modeled by adding the fluctuating (the last term) and dissipative ( γ ṙi term) forces to the Newtonian equations of motion. To integrate the Langevin equation, NAMD uses the Brünger-Brooks-Karplus (BBK) method, a natural extension of the Verlet method for the Langevin equation. The position recurrence relation of the BBK method is:

r n+1 = r n + 1 0.5 γ ∆t 1 + 0.5 γ ∆t (r n r n 1 )+ ∆t 2 1 + 0.5 γ ∆t " F(r n ) m + r 2 γ k B T m Z n # (3.17)
where Z n is a set of Gaussian random variables of zero mean and variance 1. The BBK integrator requires only one random number for each degree of freedom. The steady-state distribution generated by the BBK method has an error proportional to ∆t 2 , although the error in the time correlation function can have an error proportional to ∆t.

NPT Ensemble Simulations

Often, it is useful to maintain a simulated system at both constant temperature and pressure. In a thermodynamic sense, systems at constant pressure are the ones that can exchange volume with their surroundings (e.g., by way of a piston). Their volume therefore fluctuates. Likewise, simulated systems at constant pressure involve volume fluctuations.

A number of different barostat techniques exist with the scope of maintaining a target pressure by dynamically adjusting the volume of the system during the simulation. Mainly, the most common barostat techniques are based on:

. Volume rescaling -the instantaneous pressure is made to equal the target pressure by rescaling the system volume at periodic intervals.

. Berendsen barostat -the pressure is weakly coupled to a pressure bath and the volume periodically rescaled.

. Extended ensemble barostat (also known as Andersen barostat) -the system is coupled to a fictitious pressure bath using an extended Lagrangian and the introduction of new degrees of freedom.

A detailed discussion of the algorithms commonly used is beyond the scope of this thesis.

Initial state of the system

In order to perform a simulation for measuring some selected properties of a given system, it is essential, like in a real experiment, to prepare first the initial state of the system. This means to obtain, in some way, a set of coordinates and velocities for all the atoms that identify an initial configuration of the system compatible with the desired initial state.

Specifically, for the simulation of large biological molecules such as proteins or nucleic acids, and as a consequence of the complex structures of these macromolecules, the initial positions of their atoms are usually obtained from the results of some real experiments, like X-Rays crystallography or nuclear magnetic resonance (NMR) spectroscopy, whereas the initial velocities are typically set pseudorandomly so that the total kinetic energy of the system corresponds to the expected value at the target temperature.

Actually, the sets of coordinates obtained from these experiments often do not represent an initial configuration that is compatible with the desired initial state of the system (e.g. due to the hydration of the sample or to the fact that the set of coordinates usually refer only to one macromolecule, while the real system might be formed by more). However, since these coordinates represent the three-dimensional structure of the macromolecules involved in the simulation, they provide a fundamental starting point for the preparation of the system. Indeed from these coordinates, through several processes that exploit and modify them, it is possible to derive a configuration for the entire system that, been as close as possible to the configuration of the real system that it is intended to study, can be used to simulate its initial state.

Hence, when the initial configuration is achieved, it is possible to make several simulations of the system that can be used to measure the properties to be studied (this phase of the simulation is conventionally named: production phase and the main scope of a molecular simulation is to get some interesting results from this phase). The process generally used in MD simulations of biological system, preceding the production phase, is mainly the following: solvation, minimization, heating (set the velocities) and equilibration; that are describe below.

In any case, it is important to point out that the choice of the initial configuration must be done carefully as this can influence the quality of the entire simulation.

Solvation

Solvation consists in the process of taking into account the solvent surrounding the biomolecules. Actually, biomolecules are generally in solution with some types of aqueous solvent and therefore solvation is a common process that is used to prepare the initial configuration of the system. This is indeed due to the known fact that solvation effects play a crucial role in determining molecular conformation, electronic properties, binding energies, etc.

There are mainly two ways to solvate a biomolecular system:

• explicit treatment: the coordinates of the solvent atoms are directly added to the system with the specification of the structure of the molecules that occur in the solvent (most often water molecules and, sometimes, also salt ions).8

• implicit treatment: the force field is modified to include also the effect due to the interaction between the biomolecular system and the solvent. Thus this technique eliminates the need of specifying all the atoms of the solvent by including many of their average effects in the inter-atomic force calculation. For example, polar solvent acts as a dielectric and screens (lessens) electrostatic interactions.

The elimination of explicit water accelerates conformational explorations and sometimes increases simulation speed, although at the cost of not modeling the solvent as accurately as explicit models. However, since implicit solvent models represent the solvent in an averaged manner, they are considered less accurate than explicit solvent models. Caution should always be used when implicit solvents are employed for molecular dynamics research.

Minimization

In the context of MD simulation, the energy minimization, in short known as minimization, is the process of finding an arrangement in space for a collection of atoms where, according to the force field used, the net inter-atomic force on each atom is acceptably close to zero and the position on the potential energy surface is a stationary point. In general, the collection of atoms might be a single molecule, an ion, a condensed phase, a transition state or even a collection of any of these.

As an example, when optimizing the geometry of a water molecule, one aims to obtain the hydrogen-oxygen bond lengths and the hydrogen-oxygen-hydrogen bond angle which minimize the forces that would otherwise be pulling atoms together or pushing them apart.

The motivation for performing an energy minimization is the physical significance of the obtained structure: even when initial coordinates are available from an experiment, the starting vector may not correspond to a minimum in the potential energy function used, and as such minimization is needed to relax strained contacts. When an experimental structure is not available, a build-up technique may be used to construct a structure on the basis of the known building blocks, and minimization again is required. For a biomolecular system, the potential energy function is a very complex and multidimensional landscape. It has one deepest point, the global minimum, and a very large number of local minima. The goal of the energy minimization is to find a local minimum. The energy at this local minimum may be much higher than the energy of the global minimum. Performing an energy minimization will guarantee the removal of any strong van der Waals interactions which might otherwise lead to local structural distortion and result in an unstable simulation.

Heating

During this phase, the integration of the equations of motion starts. The heating process, usually performed in the NVT ensemble, consists of progressively increasing of the reservoir's temperature to gradually heat the system. This process is performed mainly for two reasons:

• to reach the desired temperature at which the system has to be studied;

• to change the configuration of the system so as to escape from local minima within the energy landscape that are probably obtained with the minimization but that do not represent the real state of the system. 9In both cases, the need of a gradual increase of the temperature is due to the fact that, starting from the structures obtained from the X-Ray experiment, typically performed at low temperatures, and after some manipulations and the minimization that modifies these structures, the configuration of the system obtained is built without taking into account the velocities and hence the contribution of the kinetic energy to the total energy. Indeed starting a simulation at high temperature with such artificial configuration of the system that may be considered as "frozen", might be unfeasible because the simulation will be probably numerically unstable. Therefore, during the heating phase, initial velocities are assigned at a low temperature and the simulation is started with periodically assigning new velocities at a slightly higher temperature and letting the simulation to continue. This step is repeated until the desired temperature is reached. 10

Equilibration

Once the heating process is over and the desired temperature is reached, the simulation is continued and during this phase, properties such as structure, pressure, temperature and the energy are monitored. The point of the equilibration phase is to run the simulation until these properties become stable with respect to time. If in the process, the temperature increases or decreases significantly, the velocities are scaled such that the temperature returns to its desired value.

Production and Analysis

The last step of the simulation is the production phase, wherein the system is simulated for the time length required, normally from several hundred ps to ns or more. During this process, coordinates and eventually also velocities, of the system at different times are stored in the form of trajectories. These are then used for calculations of mean energy, mean square distance (MSD) between structures, local mean square atomic fluctuations (MSF) etc. From MD simulations, time dependent properties such as correlation functions can also be calculated and these in turn can be related to spectroscopic measurements.

Coarse graining

The atomistic simulations described in the previous sections have been successfully employed to investigate many aspects of biological systems like the folding of small proteins.However, biological systems and phenomena exhibit complexities and diversities that spread over a wide and disparate range of spatio-temporal scales. These phenomena include the dynamics of large proteins and the self-assembly of biological molecules. Different space and time resolutions are involved, from the quantum mechanical level, describing the electronic structures, to the atomic scale, and the continuum level of fluid motion at macroscopic scales (Figure 3.10).

Production and Analysis

Despite the growth of computational resources and power, and the emergence of enhanced sampling methodologies, all-atom simulations are still limited to systems containing tens or hundreds of thousands of atoms on a µs time scale, that only rarely reach the millisecond time-scale [START_REF] Kmiecik | Coarse-grained protein models and their applications[END_REF]. Hence, the deployment of simplified representations still able to capture the essential features of the phenomena is required (Fig. 3.11). One successful strategy is to reduce the number of degrees of freedom by a systematic coarse-graining. Coarse-grained (CG) models present an attractive alternative to the traditional atomistic simulations offering the possibility of investigating complex cellular processes over larger lengths and longer time scales at a reduced level of detail. A coarse-graining operation requires the selection of the level of description. Then, the degrees of freedom are reduced by averaging off the behaviour of the fast ones. In the coarse-graining of biomolecules, the strategy is, for example, to cast together groups of atoms, thus preserving some degrees of molecularity. However, depending on the problem under consideration a more aggressive coarsening can be performed, i.e. a whole protein can be modelled as a single particle. The parametrization of force fields that are both accurate and transferable -that is, capable of describing the general dynamics of systems with different compositions and different configurations -becomes increasingly difficult as the graining becomes 'coarser', because more specific interactions must effectively be included in fewer parameters and functional forms. On other hand, coarse-graining also affects thermodynamic properties of a modelled system, particularly the balance between enthalpy and entropy. Reduction of the degrees of freedom affects the entropy of the simulated system, which is compensated by reduced enthalpic terms. In turn, a coarse-grained model may accurately reproduce free energy differences but contributing enthalpy and entropy values may be inaccurate. Such limitations are typical for the majority of coarse-grained models.

For fluids like water, the molecular CG consists in describing an ensemble of particles as a single entity and inserting important features at the level of the interaction potentials. Stronger simplification can reduce the fluid to essential variables, like density or velocity. For instance, in the Brownian Dynamics, the effects of collisions of the solvent with a large particle are described without representing explicitly the solvent molecules.

It is clear that when the CG is performed only phenomena at the pertinent scales are accessible and the high level information is lost. For instance, when liquid water is described at the mesoscopic level, we can not use this representation to investigate the kinetics of the hydrogen-bond formation. The reason is that the dynamics is averaged out in the model. Many strategies, and rigorous procedures can be devised to perform coarse-graining.

The OPEP force field

In our work, for the description of the proteins in the CG simulations, we have employed the Optimized Potential for Efficient protein structure Prediction (OPEP) force field. This model has been developed 21 years ago by Derreumaux and coworkers, [START_REF] Sterpone | The OPEP protein model: from single molecules, amyloid formation, crowding and hydrodynamics to DNA/RNA systems[END_REF]. and it consist in an intermediate resolution model where each amino acid is represented with six beads: the backbone is retained in full atomic detail (all N, C α , C, O and H main-chain atoms are considered), while side chains are represented by a unique bead located at the center of mass of their heavy atoms 11 , see Figure 3.12. The Van der Waals radius and the positions of the side chain beads were calculated using a database of 2,250 PDB structures with sequence identity lower than 30%. OPEP is highly transferable among different proteins, but it is limited to soluble protein in water. On the other hand, it requires a small integration time step of 1.5-2.0 fs due to the detailed backbone and hydrogen bonds, limiting the amount of sampling that can be accomplished.

Hydrodynamic interactions: Lattice Boltzmann Molecular Dynamics coupled with OPEP

Treating the hydrodynamic interactions is complex and computationally demanding, because they are long range interactions, non-linear in nature, and cannot be expressed simply as a sum of two-body terms. Several alternative schemes have been developed along the years. The basic idea is to track the solvent degrees of freedom through a simplified mesoscopic representation, with the local dynamics that satisfies the mass, momentum and energy conservation laws and recovers the solution of hydrodynamic equations in the large-scale limit. In order to simulate a soft matter system, the solvent model must then be coupled with the respective algorithms that models the dynamics of solute particles. One of the most widely used schemes is the Lattice Boltzmann (LB) method. Proteins are described at microscopic level, interact according to the OPEP coarsegrained force field and move in the continuum. The aqueous solvent is handled by the lattice Boltzmann method, so that fluid populations that reside on a Cartesian mesh are evolved in time and move to neighboring mesh points as connected by a set of discrete speeds (indicated as D3Q19) [START_REF] Chiricotto | Hydrodynamic effects on β -amyloid (16-22) peptide aggregation[END_REF].

In the LB approach, the fluid is represented through particles that reside on a threedimensional cubic lattice with spacing ∆x [START_REF] Chen | Lattice boltzmann method for fluid flows[END_REF][START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF]. Here "fluid particles" do not correspond to single molecular entities but represent instead the collective motion of the fluid (see Fig. 3.13). In short, we track the solvent kinetic degrees of freedom through a simplified representation, with a local dynamics that satisfies the mass and momentum conservation laws and recovers in the large-scale limit the Navier-Stokes equation of fluid-mechanics.

The lattice Boltzmann equation reads as follows:

f p (x + c p ∆t, t + ∆t)= f p (x, t) ω∆t • f p f eq p (x, t)+g p (x, t) (3.18)
where f p (x, t) denotes the probability distribution of finding a particle at lattice site x at time t and moving in lattice space with discrete velocity c p , ∆t is the time-step for the numerical integration. The particle displacement occurs to the first and second lattice neighbors (D3Q19) by using 18 directions plus a null one mimicking particles at rest.

The distribution f p (x, t) evolves in space and time toward the equilibrium target f eq p with the characteristic relaxation frequency ω = 1/τ. The term g p (x, t) includes the drag force F D and extra contributions F ext as the random noise encoded in the molecular dynamics. When focusing on a solute particle suspended in the fluid, this term describes essentially the particle-to-fluid back reaction. An accurate expression of g p (x,t) is given by [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice boltzmann method[END_REF]:

g p (x, t)= w p ∆t  G • c p c 2 s + (G • c p )(u • c p ) c 2 s c 2 s G • c p c 4 s (3.19)
The force G contains any external force and the exchange of momentum induced by N moving atoms, G = F D + F ext . To lowest order, it can be shown that

g p (x, t)= w p ∆t G•c p c 2 s
, whereas in practice a higher order version is needed to ensure a global second-order accuracy of the lattice Boltzmann solver [START_REF] Melchionna | Design of quasisymplectic propagators for langevin dynamics[END_REF]. Without entering into the analytical details, the method is typically extended to account for local fluctuations at the level of the stress tensor, such that fluctuating hydrodynamics are recovered.

In addition, it can be shown that local mass and momentum of the global particlefluid elements are preserved, a key condition to obtain the correct fluid dynamic behavior. It should also be remarked that the model can be easily extended to account for thermodynamic forces. The coupling between the motion of a solute particle and the fluid is based on the assumption that momenta exchange in Stokes-like fashion, thus defining a drag force between the i-th particle of mass m i and the fluid [START_REF] Ahlrichs | Simulation of a single polymer chain in solution by combining lattice boltzmann and molecular dynamics[END_REF]:

F D (R i )= m i γ i [V i ũ (R i )] (3.20)
where V i is the atom velocity, and ũ indicates the fluid velocity field distributed over the region occupied by the atom, and γ i is a friction coefficient that in principle can vary depending on the solute particle type.

In this scheme (LB coupled with the OPEP force field), the dynamics of the particles are governed by the following evolution equations for the positions R i and velocities V i of the particles:

Ṙi = V i (3.21) Vi = F C i + F D i m i + µ i (3.22)
where F C i is a conservative force describing the sum of molecular interactions as encoded in the OPEP force field. The drag force is given by eq. (3.20). As anticipated above, F D i represents the mechanical and dissipative friction exerted between a particle and the surrounding fluid. The strength of this dissipation depends on γ i , a friction parameter that can be tuned in order to alter the response time between fluid and molecular motions. Finally µ i is a white noise mimicking the effect of the thermal collisions with the molecules of the fluid, with the mean hµ i (t)i t = 0 and hµ i (t)µ i (0)i t = 2γ i k B T , being k B the Boltzmann constant and T the temperature. Equation (3.21) is integrated over the time step ∆t MD according to the symplectic position Verlet algorithm [START_REF] Melchionna | Design of quasisymplectic propagators for langevin dynamics[END_REF] and, if ∆t MD = ∆t LB , the particle and LB dynamics are updated in a synchronous way.

In conclusion, conceptually we adopted Boltzmann kinetic theory, and its numerical representation, in order to describe the solvent as a continuum in a probabilistic sense. Traditionally, in the LB formulation the lattice space ∆x which supports the fluid kinetics is defined as a representation (coarse-graining) of the collective kinetic behaviour of a group of solvent molecules. It is also accepted that in order to observe hydrodynamic behavior down to the ∆x scale, it is usually considered that the fluid mean free path should not exceed ∆x. Since in liquid water, the molecular mean free path is of the order of a few angstroms, a subnanometric lattice space can be supported in LB, allowing for the hydrodynamic behaviour to emerge at distance larger than ∆x. In this approach, the lattice grid element must not be viewed as a volumetric entity that contains a fixed number of particles, but instead a numerical support for the probabilistic description of the averaged single particle trajectories. Horbach and Succi [START_REF] Horbach | Lattice boltzmann versus molecular dynamics simulation of nanoscale hydrodynamic flows[END_REF] have shown that this strategy is effective for the simulation of nanofluids and the obtained results agree very closely with particle-based simulations.

This coupling between the LB method and the OPEP model has been implemented in MUPHY software. Temperature is a boost for cellular metabolism, but above a certain threshold, it corrupts functional processes involving proteins and causes cell death. Whether the thermal denaturation involves the whole proteome or just a subset of critical proteins in the cytoplasm is still debated. Here, we attack the problem from a preferential angle by monitoring via QENS and multi-scale simulations the dynamical state of the E.coli proteome across the cell death temperature. Above the cell death temperature, the cytoplasm experiences a dynamical slowdown caused by the unfolding of just a small number of proteins. This small fraction is sufficient to induce the gelation of the cytoplasm. From the dynamical properties, the fraction of unfolding is extracted and used to reconstruct successfully the E. coli growth rate.

Introduction

Temperature has a significant impact on cells. Notably, membranes, proteins and nucleic acids suffer in various ways from heat. The membranous integrity can be challenged resulting in the evasion of periplasmic proteins or the entrance of harmful compounds [START_REF] Russel | Lethal Effects of Heat on Bacterial Physiology and Structure[END_REF]. Proteins are mandatory for good cellular functioning, but high temperature provokes loss of conformations and denaturation. Nucleic acids are the most stable against thermal stress [START_REF] Mackey | Thermal denaturation of whole cells and cell components of Escherichia coli examined by differential scanning calorimetry[END_REF], so that their denaturation can be considered as a minor cause of cell death. Moreover, biological migrations, extinctions, genetic divergence, and speciation can all be triggered by small changes in environmental temperature [START_REF] James | Abrupt deep-sea warming, palaeoceanographic changes and benthic extinctions at the end of the Palaeocene[END_REF][START_REF] Linda C Ivany | Cooler winters as a possible cause of mass extinctions at the Eocene/Oligocene boundary[END_REF][START_REF] Andrew P Allen | Kinetic effects of temperature on rates of genetic divergence and speciation[END_REF]. A deep understanding of the cell's thermal stability is key to model the impact of climate change on microbial organism growth [START_REF] Frey | The temperature response of soil microbial efficiency and its feedback to climate[END_REF], establish theoretical boundaries for life in extreme environments [START_REF] Igor | Physics and evolution of thermophilic adaptation[END_REF] and optimize thermal based treatments for cancer [START_REF] Donald S Coffey | Hyperthermic biology and cancer therapies: a hypothesis for the "Lance Armstrong effect[END_REF]. Yet, the factors that influence the cell's thermal sensitivity are largely unknown. The proteome's thermal sensitivity has to play a key role as a determinant for most of the temperature-dependent whole-organism activities, as proteins are the most abundant and less stable biomolecules in the cell.

Different pictures have been proposed to link the degradation of the proteome to the upper limit of the cellular thermal niche, i.e. the cell's death temperature T CD .A first essential aspect is to quantify the proteome thermal stability [START_REF] Zeldovich | Protein stability imposes limits on organism complexity and speed of molecular evolution[END_REF][START_REF] Dill | Physical limits of cells and proteomes[END_REF][START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF]. On one hand a proposed theoretical model [START_REF] Ghosh | Cellular proteomes have broad distributions of protein stability[END_REF][START_REF] Dill | Physical limits of cells and proteomes[END_REF] finds that the cell death is linked to a global catastrophe of the proteome with proteins unfolding in a narrow range of temperatures near the T CD . This picture has been challenged recently by experimental investigations of E. coli lysates and cells, and based on different techniques such as limited proteolysis [START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF] or thermal proteome profile [START_REF] Mateus | Thermal Proteome Profiling in Bacteria: Probing Protein State in vivo[END_REF], combined with mass spectroscopy. According to these studies only a small set of proteins indeed unfolds at the cell death. Thermal adaptation would result from the preferential stabilization of a homologous subset of proteins, thus indicating that the heat sensitivity of cells can be explained by a small number of proteins that serve critical physiological roles.

Actually, the proteome's thermal stability is not the only physical determinant of the cell's growth rate, which is expected to depend on the rate of protein diffusion throughout the cell, the latter being often the limiting factor of the rates of cellular biochemical processes [START_REF] Paul E Schavemaker | How important is protein diffusion in prokaryotes[END_REF]. Protein diffusion depends in turn on the temperature, especially through the contribution of the intrinsic viscosity in the high-temperature range when biomolecules start to unfold. To date, the relationship between the diffusive dynamics of proteome and the thermal sensitivity of a cell has not yet been investigated, also due to the extremely difficult challenge to represent the motions of proteins in a crowded milieu cell's cytoplasm where local concentration may vary from 200 g/L up to 400 g/L [START_REF] John | Macromolecular crowding: an important but neglected aspect of the intracellular environment[END_REF]. Here, the protein diffusive dynamics is affected by several factors, such as the presence of steric barriers given by the other macromolecules, hydrodynamic and attractive interactions and spatial heterogeneity.

On these grounds, here we provide an unprecedented picture of the dynamics of the E. coli's proteome in the nanosecond time-scale, based on state-of-the-art neutron scattering spectroscopy and multi-scale molecular dynamics simulation. We show that in E. coli the global protein diffusion is a close proxy of the bacterial metabolism, with a linear Stokes-Einstein dependence in the lower temperature range and a striking dynamic slow-down above the thermal death. Combining the results on the proteome dynamics from neutron scattering and simulations we describe the way the unfolded protein fraction progressively increases with temperature, offering an alternative quantification to existing ones [START_REF] Dill | Physical limits of cells and proteomes[END_REF][START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF][START_REF] Mateus | Thermal Proteome Profiling in Bacteria: Probing Protein State in vivo[END_REF]. We clearly show that no global proteome unfolding occurs at cell death. Finally, we verify that the derived proteome stability curve and temperature-dependent proteome diffusivity together, allow to excellently reproduce the E. coli growth rate profile.

Methods

Sample Preparation

E. coli BL21 (DE3) was grown overnight in LB medium (made with H2O) at 37 C with shaking (200rpm). 3.3g of E. coli cells were collected by centrifugation and washed twice with a buffer made in deuterium oxide (99.9 atom D) as followed: cells were suspended in 36mL of D2O buffer at pD8, spun at 5400g for 18 minutes at 6 C and the supernatant was removed. The pD 8 buffer contains 50mM Tris, 150mM NaCl and 5mM KCl. To obtain a pD of 8, the pH of the buffer was adjusted to 7.6 using HCl.

QENS Experiments

We performed two experiments [START_REF] Peters | Global proteome dynamics as a proxy for cellular thermal stability[END_REF][START_REF] Di Bari | Global proteome dynamics as a proxy for cellular thermal stability[END_REF] on the cold neutron backscattering spectrometer IN16B at the Institut Laue-Langevin (ILL) [START_REF] Frick | Recent backscattering instrument developments at the ILL and SNS[END_REF], with an energy resolution of ⇡ 0.75µeV FWHM, an energy-range of |E|31µeV, and a wave-vector coverage of 0.19 Å 1  q  1.9 Å 1 . The samples were measured for 2 hours at each temperature in a temperature range from 275 K to 348 K. The data reduction (i.e. normalization to the monitor, integration over the regions of interest of the vertically position-sensitive detector tubes, calculation of the energy axis, and centering of the elastic line positions using separate Vanadium measurements) where carried out with the built-in module for IN16B of the Mantid program [START_REF] Arnold | Mantid-Data analysis and visualization package for neutron scattering and µ SR experiments[END_REF]. The subtraction of the sample holder contribution to the signals, the normalization of the detector efficiency, and the fit of the data where performed using an in-house python module that is available on github: https://github.com/DanieleDiBari/NSAnalysis. QENS Model. The fully reduced scattering function measured from the E. coli samples, S exp (q, E), can be obtained by the convolution of the theoretical scattering function S th , describing the interaction between the neutrons and the bacteria, and the instrumental resolution, R(q, E), which is determined by a vanadium sample, a completely elastic and incoherent scatterer (see section 2.3.2) [START_REF] Bée | Quasielastic neutron scattering[END_REF]:

S exp (q, E, T )=e E 2k B T • [R(q, E) ⌦ S th (q, E, T )] (4.1)
where e E 2k B T is the detailed balance factor. Since proteins represent ⇡ 55% of the bacterial dry weight (see Table 4.1) and have a percentage of hydrogen higher (50%) than any other type of macromolecules (30%) [START_REF] Jasnin | In vivo measurement of internal and global macromolecular motions in Escherichia coli[END_REF], the major contribution to S th is originating from self-diffusive dynamics of an average protein and of the bulk water present in the sample, in this case D2O (see section 2.2.3):

S th (q, E.T )=S AP (q, E, T )+φ

• S D2O (q, E, T ) (4.2)
where S AP (q, E, T ) and S D2O (q, E, T ) are, respectively, the scattering functions of the average protein and the D 2 O, and φ is a scalar factor that weights the contribution of the solvent. φ can be estimated by the product of the D 2 O amount in the sample, which constitutes about 80% of the sample mass, and the percentage of bulk water molecules in E. coli, which is almost 60%. QENS: Fitting procedure. To mitigate the problem of overfitting due to the complexity of the system and the consequent high number of unknown parameters necessary to describe the resutling QENS signal, we tried to reduce, as much as possible, the number of free parameters that can vary during the fit. To this end, we used the measurements of the PBS-D2O Buffer to fix the q-dependence of the solvent's intensity S D2O in eq. (4.2). Moreover, we employed a three step procedure for the analysis, where at each step we perform a fit of the data and we used the resulting information to improve the model used for the fit and reduce the number of free parameter.

QENS: Energy Resolution. The instrumental resolution function R(q, E) takes several parameters into account and it strongly depends on the set-up of the instrument (section 2.3.2). The red line is the best fit of R(q, E) described by the eq. ( 4.3), meanwhile the dashed lines are the different components of R(q, E). The inset shows the Full Width Half Maximum (FWHM) of R(q, E) which represents a measure of the energy resolution of the instrument.

In particular, at the IN16B spectrometer the resolution function is well described by a sum weighted by the scalar parameter F(q) of a Gaussian and a Lorentzian function [START_REF] Frick | Recent backscattering instrument developments at the ILL and SNS[END_REF]:

R(q, E)=A(q)• 8 > < > : 1 π • F(q) • γ(q) γ 2 (q)+[E E 0 (q)] 2 +(1 F(q)) • exp n [E E 0 (q)] 2 2σ 2 (q) o σ (q) p 2π 9 > = > ; +B(q) (4.
3) where |F(q)|1. B(q) reflects a background which can depend on q. Figure 4.1 shows the fit of the Vanadium data with the eq. ( 4.3) and the measured energy resolution of the instrument, ∆E FWHM , corresponding to the Full Width Half Maximum (FWHM) of eq. ( 4.3). The resulting parameters averaged over q are reported in table 4.2. hF(q)ih A(q)ih E 0 (q)ih B(q)ih σ (q)ih γ(q)ih FWHMi QENS: D2O Contribution. To simplify the fit of the E. coli data, we first analyzed the QENS data of the PBS-D2O buffer at three temperatures: 276 K, 300 K, and 324 K to obtain an estimation of the q-dependence of the intensity of the signal due to the solvent. The main contribution here comes from D2O bulk water whose scattering function can be well represented as follows [START_REF] Grimaldo | Highresolution neutron spectroscopy on protein solution samples[END_REF]:

S D2O (q, E)=e E 2k B T • R(q, E) ⌦ ⇥ β D2O (q) • L γ D2O (E) ⇤ (4.4)
where β D 2 O is the intensity of the signal, and L γ D2O is a Lorentzian function that de- scribes the diffusive translational motions of the D2O molecules. An example of the resulting parameters is reported in Fig. 4.2. The D2O cross section comprises 79% coherent and 21% incoherent scattering and, applying Vineyard's convolution approximation, we have S D2O (q, E = 0) ⇡ S (coh) D2O (q), where S (coh) D2O is the coherent static structure factor [START_REF] George | Scattering of slow neutrons by a liquid[END_REF]. Hence, S D2O (q, E) has as main contribution the coherent scattering of D2O that can not be avoided [START_REF] Wanderlingh | Multi-component modeling of quasielastic neutron scattering from phospholipid membranes[END_REF]. It is important to observe that the intensity of S D2O (q, E) that we have measured (Fig. 4.2) reproduces qualitatively well the shape of S (coh) D2O (q) measured by Bosio et al. [START_REF] Bosio | Enhanced density fluctuations in water analyzed by neutron scattering[END_REF].

QENS: Average Protein Contribution. The dynamic scattering function of the proteins is modeled here by the following expression (see section 2.2.3) [START_REF] Jasnin | In vivo measurement of internal and global macromolecular motions in Escherichia coli[END_REF][START_REF] Grimaldo | Hierarchical molecular dynamics of bovine serum albumin in concentrated aqueous solution below and above thermal denaturation[END_REF]:

S AP (q, E, T )=I(q, T )•{L(E; γ G (q, T )) ⌦ [A 0 (q, T ) δ (E)+(1 A 0 (q, T )) L(E; γ L (q, T ))]}
(4.5) where I(q) is the intensity of the signal which is related to the vibrational modes of the proteins, meanwhile L(E; γ G (q, T )) and L(E; γ L (q, T )) are two Lorentzian functions accounting for the diffusive contributions due to the global motions (roto-translations of the entire proteins) and local dynamics (internal motions of sub-parts of the proteins, e.g. conformational changes). Finally, A 0 (q), which multiplies the delta function δ (E), represents the EISF containing information on the geometry of the internal motions that are confined with respect to the space-time window of the spectrometer [START_REF] Bée | Quasielastic neutron scattering[END_REF].

QENS: Data analysis. As described in the previous paragraphs, we employed a three step-procedure for the fit of the data.

In the 1 st step, we performed a simple fit that treated the measured spectra independently, apart for the weighting factor for the solvent contribution φ which was shared among the spectra since it should depend only on the amount of bulk D2O in the samples. On the contrary, all the remaining parameters can vary both with q ant T . As shown in Fig. 4.3 and 4.4, this allows us to study the trends of the Lorentzian widths as function of q 2 , and we verify that, both for the global and the local motions, they are well described by a jump-diffusion model as already observed in previous in vivo QENS experiments on bacteria (see section 2.2.3) [START_REF] Jasnin | In vivo measurement of internal and global macromolecular motions in Escherichia coli[END_REF][START_REF] Martinez | High protein flexibility and reduced hydration water dynamics are key pressure adaptive strategies in prokaryotes[END_REF][START_REF] Mamontov | Microscopic diffusion processes measured in living planarians[END_REF]:

γ i (q, T )= q 2 D i (T ) 1 + q 2 D i (T ) τ i (T ) for i 2{G, L} (4.6)
Within this model, the diffusion is assumed to occur via infinitely small, elementary jumps characterized by a negligible jump time during which the particle diffuses and the residence time τ i , i.e. the time a proton spends in a given position. In the 2 nd step, we performed a simultaneous fit of the spectra measured at different q-values assuming a-priori the jump-diffusion model for both the global and the local dynamics, i.e. constraining the q-dependence of γ G and γ L with the eq. (4.6). The resulting values of the diffusion coefficients and the residence time for the global and the local dynamics are reported in Fig. 4.5.

D G (T )=(a 1 T + b 1 ) h 1 a QENS u (T ) i +(a 2 T + b 2 ) a QENS u (T ) (4.7)
with:

a QENS u (T )= 1 1 e (T T 0 )/∆T (4.8)
where a QENS u is a smeared step function describing the transition between the initial liquid state and the final gel-like state.

In the 3rd and last step, we constrained the temperature dependence of D G with eq. (4.7), and we performed a simultaneous fit of all the measured spectra. The resulting parameters constant in T and q are reported in table 4.3.

Table 4.3:

Parameters shared among all the spectra (i.e. constant in T and q) obtained from the simultaneous fit of the E. coli data at all the temperatures with eq. (4.1), taking into account eq. (4.2), eq. (4.4), eq. (4.5), and assuming a priori the models described by eq. (4.6) and eq. (4.7). As it was anticipated in section 2.2.3, the EISF can be described quite well by the model proposed by Grimaldo et al. in previous study on g-globulins [START_REF] Grimaldo | Diffusion and dynamics of γ-globulin in crowded aqueous solutions[END_REF]:

A 0 = A 0 (q, T )=p L (T )+[1 p L (T )] ⇥ s L (T ) A sph (q, T )+(1 s L (T )) A 3JD (q, T ) ⇤ (4.9)
where:

A sph (q, T )= ✓ 3 • j 1 (q • r(T )) q • r L (T ) ◆ 2 and A 3JD (q, T )= 1 + 2 j 0 (q • a M ) 3 (4.10)
p L is the fraction of H-atoms appearing fixed on the accessible time scale of the instrument. For the remaining H-atoms, two types of motions were taken into consideration: confined diffusion in an impermeable spherical volume of radius r L , described by the amplitude A sph , and random jump diffusion between three equidistant sites on a circle of radius a M , arising from the methyl-groups, and described by A 3JDa M is fixed to 1.715 Å which is the average distance of H-atoms in methyl groups. The relative contribution of these two diffusive processes is measured by s L , which, more specifically, describes the fraction of H-atoms that, among those atoms that are not appearing fixed (i.e. 1 p L ), are undergoing the spherically confined diffusion.

EINS Experiments

As Elastic Incoherent Neutron Scattering (EINS) measurements are much faster, we performed also EINS measurements for much more temperature points on the E. coli samples. The elastic data were collected for 30 seconds every 6.5 minutes, meanwhile the sample was first heated from 300K to 350K with a heating rate of 0.25K/min and then cooled down, in 20 minutes, from 350K to 305K [START_REF] Di Bari | Global proteome dynamics as a proxy for cellular thermal stability[END_REF]. In general, incoherent neutron scattering on a biological sample probes the selfdiffusion, therefore the dynamics of hydrogen atoms, which present the highest cross section among the particles contained in living matter [START_REF] Varley | Neutron scattering lengths and cross sections[END_REF]. As they are mainly distributed homogeneously in the sample, they allow to probe average molecular dynamics of the H nuclei and of the molecular groups, to which they are bound. Elastic scattering refers to processes within the time limit going to infinity, therefore only short local motions can be resolved. The elastic signal corresponds to the static structure function S(q, ∆E) , which can be written approximately as polynomial in q 2 [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF] containing the atomic mean square displacement (MSD) ⌦ u 2 ↵ :

S(q, E ⇡ 0)=e 1 6 •[a + hu 2 iq 2 + cq 4 ] (4.11)
Taking the logarithm of this expression gives access to ⌦ u 2 ↵ as function of temper-ature, which can then be fitted by applying the same gelation model as the one used in eq. (4.7) to fit the global diffusion coefficient. Fig. 4.7 shows the resulting MSD and their fit.

Model preparation for the simulations

Protein composition. To represent the protein composition of the E. coli cytoplasm (see section 3.2.4), we built our simulation system on the basis of a previous computational model [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], which was derived from the results of a proteomics study [START_REF] Andrew | Comparing the predicted and observed properties of proteins encoded in the genome of Escherichia coli K-12[END_REF] of E. coli grown under minimal media conditions. The computational model reported in [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF] contained 45 different protein species as well as 5 types of RNA and RNA-protein complexes. For the purposes of the present study and in order to permit back-mapping of smaller sub-volumes into an all-atomistic resolution, we reduced the model in the following aspects. First of all, we considered a smaller cytoplasmic volume than [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], namely a 400 Å cubic box. Second, we focused exclusively on proteins since they are known to form the most abundant macromolecular type in the E. coli cytoplasm and since the second most abundant macromolecular type -the RNA -is mainly concentrated in large ribosomal particles, which we did not include in our model. Third, we only considered protein species with molecular weights below 150 kDa to avoid large protein oligomers that would not fit easily in back-mapped sub-volumes. We assumed the same target macromolecular concentration as [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], that is, 275 g/L. Analogously to [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], the number of copies of each individual protein species was based on the relative abundances reported in [START_REF] Andrew | Comparing the predicted and observed properties of proteins encoded in the genome of Escherichia coli K-12[END_REF]. From the list of proteins simulated by [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], we omitted those counting less than 0.5 copies per our simulation box. In addition, we did not include the GltD, Hns, Pnp, and Bcp proteins, each of which would only exist in one or two copies in our simulation box, owing to the absence of reliable structures for homology modeling. On the other hand, for computational investigations of protein unfolding that are not reported in this work, we included 10 small monomeric barrels of superoxide dismutase 1 (SOD1) [START_REF] Danielsson | Global structural motions from the strain of a single hydrogen bond[END_REF], raising the total macromolecular concentration to 279 g/L. Overall, the system comprised 197 proteins of 35 species (see Table 4.4 and Figure 4.8).

Obtaining protein structures. Where available, we used an E. coli PDB structure to prepare an atomistic model of each protein species (see section 3.2.4). Unless we found a more recent and higher-quality structure in the PDB database [START_REF] Berman | Announcing the worldwide protein data bank[END_REF], we selected the same PDB structure as was used in the previous work [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF]. In four cases, we made use of a homology model stored in the SWISS-MODEL repository [START_REF] Waterhouse | SWISS-MODEL: homology modelling of protein structures and complexes[END_REF] as a starting point of the model preparation. Where needed, we subsequently added missing residues and corrected non-proline cis peptide bonds using the MODELLER software, version 9.23 [START_REF] Šali | Comparative protein modelling by satisfaction of spatial restraints[END_REF]. The quality of the structure was checked with the MolProbity server [START_REF] Christopher | MolProbity: More and better reference data for improved all-atom structure validation[END_REF].

When present in the PDB structure, metal ions coordinated to the protein were included in the model; however, larger ligand molecules were omitted. With the exception of TufA and GpmA, the oligomerization state of the protein was taken to be the same as in [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF]. While TufA was modeled as a dimer by [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], the prediction of the PDBe PISA server [START_REF] Krissinel | Inference of macromolecular assemblies from crystalline state[END_REF] yielded an ambiguous result, and the protein was described as a monomer in a previous experimental work [START_REF] Wiborg | Mapping Escherichia coli elongation factor Tu residues involved in binding of aminoacyl-tRNA[END_REF]. Therefore, we modeled TufA in a monomeric state. Similarly, while the GpmA protein was simulated in a dimeric form by [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF], the PDBe PISA prediction placed it in a grey zone, and according to the EcoCyc database [START_REF] Keseler | EcoCyc: a comprehensive database of Escherichia coli biology[END_REF], the ATP-free version of the protein was monomeric. As a result, we considered GpmA in a monomer state. The atomistic protein structures were hydrogenated by the GROMACS pdb2gmx tool [START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF], placed in rectangular boxes with a minimum distance of 1.5 nm between the protein and the boundary of the box, and solvated in a 150 mM KCl solution, including several additional ions to neutralize the net charge of the system, if needed. Each system was then subjected to a short energy minimization and equilibration protocol, performed in the GROMACS 2018.7 software [START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF]. First, an energy minimization was performed to bring the maximum force below 1000 kJ mol 1 nm 1 , while harmonic restraints were applied to all heavy atoms of the protein (with a force constant k BB = 1000 kJ mol 1 nm 2 assigned to the backbone atoms and a force constant k SC = 500 kJ mol 1 nm 2 applied to the side-chain atoms). This energy minimization was followed by a six-step relaxation protocol with gradually decreasing harmonic position restraints, where the first two short simulations were performed in the NV T ensemble and were followed by four NPT trajectories simulated at a pressure of 1.01 bar (see Table 4.5 for more details). The minimization and equilibration protocol was performed separately for the Amber99SB-disp [START_REF] Robustelli | Developing a molecular dynamics force field for both folded and disordered protein states[END_REF] and CHARMM36m [START_REF] Huang | CHARMM36m: an improved force field for folded and intrinsically disordered proteins[END_REF] force fields. The same simulation parameters were used as for the respective Amber99SB-disp and CHARMM36m systems described in the subsection 4.2.7, except for the temperature coupling, which was ensured by the Berendsen thermostat [START_REF] Herman | Molecular dynamics with coupling to an external bath[END_REF] with a time constant of 1 ps. The equilibrated protein structures obtained with the Amber99SB-disp force field (Figure 4

LBMD simulation

We performed a coarse-grained simulation of the large system using the Lattice Boltzmann Molecular Dynamics (LBMD) approach (see section 3.3.3) [START_REF] Sterpone | Protein simulations in fluids: Coupling the OPEP coarse-grained force field with hydrodynamics[END_REF], implemented in-house in the MUPHY software [START_REF] Bernaschi | MUPHY: A parallel MUlti PHYsics/scale code for high performance bio-fluidic simulations[END_REF]. Our LBMD scheme, which has been successfully applied to a number of biological systems [START_REF] Chiricotto | Hydrodynamic effects on β -amyloid (16-22) peptide aggregation[END_REF][START_REF] Chiricotto | Multiscale aggregation of the amyloid Aβ 16-22 peptide: from disordered coagulation and lateral branching to Amorphous prefibrils[END_REF][START_REF] Sterpone | Molecular mechanism of protein unfolding under shear: A lattice boltzmann molecular dynamics study[END_REF][START_REF] Brandner | Modelling lipid systems in fluid with Lattice Boltzmann Molecular Dynamics simulations and hydrodynamics[END_REF], combines a coarse-grained protein model with a lattice-based description of hydrodynamic interactions [START_REF] Sterpone | Protein simulations in fluids: Coupling the OPEP coarse-grained force field with hydrodynamics[END_REF]. The protein model was based on the OPEP v.4 force field [START_REF] Sterpone | The OPEP protein model: from single molecules, amyloid formation, crowding and hydrodynamics to DNA/RNA systems[END_REF]. Since the goal of this coarse-grained simulation was to sample different intermolecular arrangements rather than simulate protein conformational changes, the conformations of the proteins were restrained by an elastic network (distance cutoff 6 Å, force constant 5 kcal mol 1 Å 2 ). The primary focus on intermolecular interactions also allowed us to simplify the description of the protein backbone, which was exclusively represented by C al pha beads, the sizes of which were increased by 50 % relative to the standard OPEP v.4 model. Finally, the OPEP sidechain-sidechain non-bonded interactions were rescaled by a factor of 0.857 to mitigate excessive aggregation [START_REF] Timr | Stabilizing or Destabilizing: Simulations of Chymotrypsin Inhibitor 2 under Crowding Reveal Existence of a Crossover Temperature[END_REF].

The simulation was performed in an NV T ensemble at T = 300 K, using a time step of 10 fs for bonded interactions and a time step of 20 fs for non-bonded and hydrodynamic interactions. A trajectory length of 4.3 µs was reached. Hydrodynamic interactions were described using the LB technique [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF], employing the BGK (Bhatnagar-Gross-Krook) collisional operator [START_REF] Benzi | The lattice Boltzmann equation: theory and applications[END_REF], with a lattice grid spacing of 4 Å. The LB kinematic viscosity ν 0 was set to reproduce bulk water behaviour at ambient condi- tions. The coupling coefficient γ was determined for each protein species as a function of its molecular weight (see Figure 4.9). The dependence of γ on molecular weight was derived by performing separate simulations of six protein species (namely CspC, SOD1, Adk, TufA, MetE, and GapA) in dilute conditions and by adjusting the γ of each species to match the diffusion coefficient of each isolated species with a prediction obtained with the HYDROPRO software [START_REF] Ortega | Prediction of hydrodynamic and other solution properties of rigid proteins from atomic-and residuelevel models[END_REF].

Sub-volume selection and back-mapping

To explore the diffusion of proteins with an all-atom resolution, allowing the description of conformational flexibility and unfolding, we selected five sub-volumes sampled in the LBMD trajectory and converted them into an atomistic resolution.

The selection of the sub-volumes proceeded in the following way. A cube with an edge length of 170 Å was placed in randomly chosen trajectory frames over the first 3 µs of the LBMD trajectory, with random positions and orientations (60,000 trials). In the dense environment of the cytoplasmic system, it was practically impossible to find sub-volume positions and orientations so as to avoid cutting at least a few proteins by the sub-volume boundary. For each placement, proteins that were fully contained in the sub-volume were counted as well as those that were cut by the sub-volume boundary. Subsequently, we isolated 120 sub-volumes that minimized the ratio of the total mass of proteins cut by the sub-volume boundary versus the total mass of proteins entirely placed inside the sub-volume. If we had only kept those proteins that were not cut by the sub-volume boundary, we would have systematically underestimated the protein concentration in the sub-volumes. Therefore, for each of the 120 sub-volumes, the proteins that were integrally inside the sub-volume were complemented by a subset of proteins that were cut by the boundary. This subset was constructed in the following way. First of all, we only retained proteins the periodic images of which did not overlap with any protein that was integrally placed inside the sub-volume. An overlap was defined as a close contact between two protein beads, with a distance less than 3 Å. Next, we checked for overlaps among the proteins that we retained in the previous step, and we removed, one by one, proteins that exhibited the largest number of overlaps with the remaining members of the subset. In this way, we ended up with a subset of proteins that had no overlap with any other protein in the sub-volume.

Out of the 120 sub-volumes constructed with the procedure described above, we selected a representative set of five sub-volumes that 1) approximated the distribution of local protein concentrations inside the large volume and that 2) maximized the diversity of protein species represented in the set while approximating the number distribution of protein species in the large volume. The distribution of local protein concentrations (see Figure 4.21B) was obtained by randomly placing a 170 Å cube over the trajectory (2000 trials) and by evaluating the protein concentration inside the cube for each placement.

The protein concentration and species composition of each of the five cubic subvolumes, sampling the structure and local concentration of the crowded solution, is indicated in Table 4.6. Each sub-box, with a size of 170 Å and containing between 11 and 27 proteins, was subsequently converted into the all-atom resolution using the following back-mapping protocol. Atomistic structures prepared as described in the paragraph "Obtaining protein structures" were overlapped with the geometries from LBMD by aligning the C al pha atoms of the atomistic structures with the corresponding C al pha beads. We checked the back-mapped configuration for the presence of entangled aromatic side chains and corrected this artifact if needed. The boxes were subsequently hydrated, and a 150 mM concentration of K + and Cl ions was added, including extra ions neutralizing the net charge of the proteins. The systems were then minimized and equilibrated using the protocol described below.

All-atom simulations

Force field and simulation parameters. All-atom molecular dynamics (MD) simulations were performed using the GROMACS 2019.4 software [START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF]. We employed two distinct sets of force field parameters to describe the proteins (see section 3.2.1): Amber Amber99SB-disp [START_REF] Robustelli | Developing a molecular dynamics force field for both folded and disordered protein states[END_REF] and CHARMM36m [START_REF] Huang | CHARMM36m: an improved force field for folded and intrinsically disordered proteins[END_REF], which belong to the most recent generation of protein force fields and which were optimized to correctly capture the properties of unfolded ensembles. The Amber99SB-disp protein force field was coupled with the Amber99SB-disp water model, while CHARMM36m used the TIP3P water model [START_REF] William L Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF]. In both cases, K + and Cl ions were described with the default parameters for the respective force field. A 1.2 nm cutoff was applied to short-range non-bonded interactions. In addition, van der Waals forces were smoothly switched to zero between 1.0 and 1.2 nm in the CHARMM36m simulations. Long-range electrostatic interactions were evaluated using the particle mesh Ewald method [START_REF] Darden | Particle mesh Ewald: An N log(N) method for Ewald sums in large systems[END_REF]. Newton's equations of motion were propagated using the leap-frog algorithm [START_REF] Williams Hockney | Quiet high-resolution computer models of a plasma[END_REF]. The LINCS algorithm [START_REF] Hess | LINCS: a linear constraint solver for molecular simulations[END_REF] was used to constrain the lengths of all protein bonds involving hydrogen, and the SETTLE algorithm [START_REF] Miyamoto | Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models[END_REF] was employed to keep water molecules rigid.

The temperature of the system was maintained by the velocity rescaling thermostat with a stochastic term [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF], which was coupled separately to the proteins and to the rest of the system with a time constant of 0.1 ps. The pressure was kept at 1.01 bar by the Parrinello-Rahman barostat [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] in all the production runs and by the Berendsen barostat [START_REF] Herman | Molecular dynamics with coupling to an external bath[END_REF] in all the NPT equilibration steps, with the time constant equaling 2 ps for both barostats.

Sub-volume equilibration. The procedure that we employed to equilibrate the atomistic sub-volumes at different temperatures and both in folded or unfolded states is detailed in Figure 4.10). First, the configurations back-mapped from LBMD underwent a short initial relaxation (Figure 4.10B). Subsequently, a heating protocol (Figure 4.10C) was used to equilibrate folded sub-volumes at increased temperatures, ranging between 310 and 380 K. In parallel to this, an unfolding protocol (Figure 4.10D) and a subsequent cooling protocol (Figure 4.10E) were employed to produce unfolded sub-volumes and equilibrate them at different temperatures in the range between 300 and 380 K. To efficiently heat up or cool down the system in the heating and cooling protocols, a heating sequence was used [START_REF] Gavin | Accelerated simulation of unfolding and refolding of a large single chain globular protein[END_REF], enabling us to quickly reach the desired temperature while allowing the system to equilibrate. To follow the unfolding of the proteins during the unfolding protocol, we measured the root-mean-square deviation (RMSD) of each protein's atoms, its radius of gyration (R g ), as well as the amount of secondary structure present in the system, as determined by the DSSP software (Figure 4.11) [START_REF] Wouter G Touw | A series of PDB-related databanks for everyday needs[END_REF][START_REF] Kabsch | Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features[END_REF].

Preparation of partially unfolded sub-volumes. The partially unfolded sub-volumes were derived from the same initial configuration as the 287 g/L sub-volume (Table 4.6). The proteins selected for unfolding are listed in Table 4.7. The unfolded fraction r u was quantified by computing the fraction of hydrogen atoms belong to the unfolded protein structures. The sub-volume preparation followed a protocol equivalent to that employed to obtain fully unfolded sub-volumes (Figure 4.10). In addition, the positions of the atoms in proteins deemed to remain folded were kept frozen during the unfolding protocol (Figure 4.10D) as well as during the three pre-cooling stages of the cooling protocol (Figure 4.10E). Moreover, before the subsequent 100 ns NPT equilibration at 360 K, the frozen folded structures were gradually relaxed using the same six-step NV T equilibration as in the initial relaxation protocol (Figure 4.10B); however, in contrast to the initial relaxation, the decreasing harmonic position restraints acted on the atoms of the folded structures only.

Production simulations. The production runs following the equilibration procedure (see section 3.3.1) were performed at eight different temperatures (300, 310, 320, 330, 340, 350, 360, and 380 K) both for folded and unfolded sub-volumes with a time step of 2 fs and each reaching a trajectory length of 102.4 ns. For the 288 g/L sub-volume simulated using the CHARMM36m force field, the production trajectories performed at T = 330 K were extended to 1 µs for further analysis (see Figure 4

Quick energy minimization with position restraints

INITIAL RELAXATION

End when Fmax < 1000 kJ mol -1 nm -1

Force constants of harmonic position restraints: kbb = 2500 kJ mol -1 nm -2 (backbone heavy atoms) ksc = 1000 kJ mol -1 nm -2 (side chain heavy atoms) Obtaining average D t per sub-volume. To evaluate the average protein translational diffusion coefficients D t for a given sub-volume, folding state, and temperature, the production trajectory was divided into 20 ns blocks. In each block, we first calculated the mean squared displacement (MSD) of the center of mass of each protein molecule.

Back-mapped configuration

Relaxed configuration

We then computed the average of these MSD curves weighted by the numbers of atoms of the proteins (see Figure 4.12). This average MSD curve was subsequently fitted with a straight line in the 0.3-5 ns regime, and D t was extracted from the slope of the fit. Subsequently, we determined the average of D t over the 20 ns blocks as well as the corresponding standard error of the mean. For the given sub-volume, folding state, and temperature, the resulting D t was corrected for the effects of periodic boundary conditions (PBC) (see the section 4.2.9). The estimated error of the PBC correction was then added to the standard error of the mean determined from the block averaging to express the uncertainty of D t .

To obtain a final value of D t characterizing the entire cytoplasm, we calculated an average of the PBC-corrected D t 's over all the sub-volumes, with the individual D t 's being weighted by the numbers of protein hydrogen atoms in the sub-volumes. An analogous weighted average of the error bars was performed to estimate the uncertainty of the final value of D t . Calculating D t and D r for the evaluation of D app . To compute the apparent diffusion coefficient D app (see section 2.2.3), a quantity directly comparable with the experimental observable, it was necessary to calculate D t and the rotational diffusion coefficient D r separately for each individual protein chain. The reason for considering individual chains rather than entire protein molecules lies in the fact that upon unfolding, the different sub-units of a protein may dissociate, and as a consequence, the evaluation of D r for the entire molecule would no longer be meaningful.

For each protein chain, D r was computed by using the following approach [START_REF] Wong | Evaluating rotational diffusion from protein MD simulations[END_REF][START_REF] Nawrocki | Slow-Down in Diffusion in Crowded Protein Solutions Correlates with Transient Cluster Formation[END_REF], approximating the chain rotation as being isotropic. One thousand randomly distributed unit vectors were centered and rotated along with the protein chain. This was done -for each frame -by fitting the geometry of the protein chain to a reference structure to which the fixed orientations of the unit vectors had been set. An autocorrelation function of the unit vectors' directions was then calculated by employing the gmx rotacf tool [START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF] and using the second-order Legendre polynomial for evaluating the autocorrelation. The resulting autocorrelation curve was fitted with an exponential function in the 0.3-5 ns regime to obtain the decay time τ which was subsequently used for calculating D r as D r = 1/(6τ).

To achieve the best possible convergence, the single-chain D t and D r were calculated from the entire trajectory without cutting it into multiple blocks.

Correcting D t and D r for PBC effects

To correct the translational diffusion coefficients for PBC effects, we added the following term [START_REF] Yeh | System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions[END_REF] to D t computed from simulations:

D corr t = k B T ξ 6π • η(T, c) • L (4.12)
where ξ = 2.837297, L is the side length of the cubic box, and η(T, c) is an empirical function that we parameterized on the basis of viscosity values computed from simulations (this procedure is described in details in the following section). The uncertainty in the value of η(T, c) was propagated to estimate the error of the correction term D corr t . Similarly, we computed the correction to the rotational diffusion coefficient as [START_REF] Linke | Rotational Diffusion Depends on Box Size in Molecular Dynamics Simulations[END_REF]

D corr r = k B T 6 • η(T, c) • L 3 (4.13)

Viscosity calculations

We estimated the low-frequency, low-shear viscosity η of the 288 g/L sub-volume (see Table 4.6) at three different temperatures (300 K, 320 K, and 380 K) for both the folded-and the unfolded systems. We also considered the intermediate unfolded systems with 25%, 50%, and 75% of unfolded proteins (see Table 4.7) at T = 300 K. All the calculations were repeated for both sets of force field parameters employed in this study.

To calculate η, we followed a similar approach as was used in [START_REF] Sören Von Bülow | Dynamic cluster formation determines viscosity and diffusion in dense protein solutions[END_REF]. For each system and temperature, we extracted 30 snapshots (50 for the unfolded CHARMM36m system at 300 K) in 1 ns intervals from the final part of the respective production run. Each of those snapshots served as a starting point of a 1 ns NPT equilibration followed by a 10-30 ns NV T simulation (10 ns for the folded systems, 20 ns for the rest of the Amber99SB-disp systems, and 30 ns for the remaining CHARMM36m systems). The NV T simulation was used to sample the elements of the pressure tensor, which was saved in 10 fs intervals. The viscosity was calculated from an integral of the autocorrelation functions of the pressure tensor fluctuations. Namely, for each off-diagonal element of the pressure tensor (P ij = P xy , P xz , and P yz ) and for each of the three combinations P ij =(P xx P yy )/2, (P xx P zz )/2, and (P yy P zz )/2 of its diagonal elements, we computed the running integral

η ij (t)= V k B T Z t 0 hP ij (0)P ij (τ)i dτ (4.14)
where V represents the volume of the simulation box. For each ij, the average of η ij (t) over all the NV T simulations was fitted with an analytical function corresponding to a bi-exponential decay of the autocorrelation function:

η ij (t)=A ατ 1 (1 e t/τ 1 )+B (1 α) τ 2 (1 e t/τ 2 ) (4.15)
where A, B, α, τ 1 , and τ 2 were parameters of the fit. The cutoff time t cut for the fitting was determined as the time where the standard error of the mean of η ij (t) for the given ij exceeded 20% of the mean. Examples of η ij (t) and the fits can be found in figure 4.13. The η ij was then calculated as the limit of η ij (t) at infinity. Finally, the value of η was obtained as the mean of the six η ij values, and the uncertainty of η was determined as the standard error of the mean of the six η ij values.

To obtain a dilute reference, we also calculated the viscosity of an aqueous solution of 150 mM KCl considering the two respective water models-the Amber99SB-disp water model and the TIP3P model-and the same three temperatures as above. In each case, we performed a 100 ns NV T simulation of a 5.1 nm cubic box after a short energy minimization and a 1 ns NPT equilibration. We split the simulation into 10 ns blocks and computed the autocorrelation functions hP ij (0)P ij (t)i separately for each block. Next, we averaged the autocorrelation functions over ij for each block and calculated the running integral η(t) from the average. Subsequently, the average of η(t) from all the blocks was fitted-using a 10 ps cutoff-with the same analytical function as shown in Eq. (4.15). Analogously to the case of the crowded sub-volumes, η was computed as the limit of the analytical function at infinity. The η ij (t) running integrals (Eq. (4.14)) obtained for the 288 g/L CHARMM36m sub-volume at T = 300 K. The black curve shows the average of η ij (t) over all the NV T simulations, the gray stripe represents the standard error of the mean, and the blue curve is a fit using the analytical function (4.15).

To estimate the correction to the diffusion coefficients due to periodic boundary conditions, we interpolated the viscosity to intermediate temperatures and extrapolated it to different protein concentrations by using analytical expressions derived in [START_REF] Monkos | Viscosity of bovine serum albumin aqueous solutions as a function of temperature and concentration[END_REF]. First, we interpolated the temperature dependence of the solvent viscosity η 0 (T ) using the following expression [START_REF] Monkos | Viscosity of bovine serum albumin aqueous solutions as a function of temperature and concentration[END_REF]:

η 0 (T )=exp ✓ B w + D w T + ∆E w RT ◆ (4.16)
The curves resulting from these interpolations are shown in As a next step, we interpolated the viscosities of the crowded system to different temperatures using the relationship [START_REF] Monkos | Viscosity of bovine serum albumin aqueous solutions as a function of temperature and concentration[END_REF]:

η(T, c)=η 0 (T ) exp  c α β c ✓ B + DT + ∆E RT ◆ (4.17)
In this expression, c denotes the protein concentration while B, D, and ∆E are free parameters. The parameter α is a function of the protein molecular weight M p , the molecular weight of water M w , and its density ρ w :

α = ρ w M p M w (4.18)
We set M p to the average molecular weight of the proteins in our model. The parameter β was defined as:

β = αν 1 (4.19)
where we set ν to a value that was determined previously for BSA, ν = 1.417 x 10 3 m 3 /kg [START_REF] Monkos | Viscosity of bovine serum albumin aqueous solutions as a function of temperature and concentration[END_REF] The parameters B, D, and ∆E obtained for each force field and for both folded and unfolded systems are listed in Table 4.9 and the resulting η(T, c) functions for c = 288 g/L are displayed in Fig. 4.15A, and Fig. 4.15B. As illustrated in Figure 4.15C, the analytical function η(T, c) defined by Eq. (4.17) allowed us to extrapolate the dependence of viscosity on concentration for a given T . To estimate the uncertainty of the interpolated/extrapolated viscosity values, we repeated the fitting of Eq. (4.17) with random combinations of viscosities each sampled within the confidence interval at the given temperature.

Evaluation of the apparent diffusion coefficient

In order to directly compare simulations with experiments, we estimated the apparent diffusion coefficient D of each protein, which takes into account not only the translational motions, but also the rotations of the entire protein. To this end, we used the following relation that links D to the translational, D t , and the rotational. D r diffusion coefficients [START_REF] Roosen-Runge | Protein selfdiffusion in crowded solutions[END_REF][START_REF] Grimaldo | Highresolution neutron spectroscopy on protein solution samples[END_REF]:

N ∑ n=0 B n (q) • D r n(n + 1)+(D t D) q 2 [D r n(n + 1)+(D t + D) q 2 ] 2 ! N!∞ 0 (4.20)
where B n (q) is:

B n (q)=(2n + 1) Z ∞ 0 ρ H (r) • j 2 n (qr) dr (4.21)
and j n is the spherical Bessel function of order n and ρ H is the Radial Distribution Function (RDF) of the protein H atoms. After the calculation of D t and D r for each protein (as described in the previous subsections), we fixed N = 550 [START_REF] Grimaldo | Highresolution neutron spectroscopy on protein solution samples[END_REF] and we solved numerically the l'eq. (4.20) for three values of q:0 .2 Å 1 ,1 Å 1 , and 2 Å 1 , verifying that, in the q-range explored by the experiment (0.2 Å 1  q  2 Å 1 ), the value of D is constant. Indeed, if N and q are big enough, it is possible to show that the result do not depend on them [START_REF] Grimaldo | Highresolution neutron spectroscopy on protein solution samples[END_REF]. In this way, we obtained the D of each protein in the folded and the unfolded configuration at all the temperatures. Then, aiming to compare these data with the QENS results, for each system we calculated the average value of D weighting each protein with the number of its H atoms. Figure 4.16 shows the results for the 5 selected subvolumes extracted from the larger cubic box used for the LBMD simulations. To combine the information related to the systems with fully folded proteins with the ones with fully unfolded proteins, we used the following relation:

D(T )=D ( f ) • [1 a QENS u (T )] + D (u) • a QENS u (4.22)
for further details on the validity of this approach see the section 4.4.2. The resulting values of D, for CHARMM36m, are reported in Fig. 4.17.

Finally, in order to calculate the most representative value of D app for the average protein in the E. coli cytoplasm, we considered the probability distribution p(c) for the selection of a sub-volume of (17 nm) 3 and concentration c from the larger system of (40 nm) 3 used in the LBMD simulation as representative of the bacterial cytoplasmsee Fig. 4.18A. The weights of each sub-volume (Fig. 4.18B) was calculated as follow:

w(c i )= p(c i ) ∑ 5 k=1 p(c k ) (4.23)
The result, obtained with both the two force field, is shown in Figure 4.28.

Results

Characterization of the dynamical state of the E. coli cytoplasm near cell death.

Quasi-elastic Neutron Scattering (QENS) experiments were performed on the backscattering spectrometer IN16b [START_REF] Frick | Recent backscattering instrument developments at the ILL and SNS[END_REF] at the Institut Laue Langevin (ILL), France, on in vivo E. coli samples. The scattering signal of these samples is mainly due to the large incoherent scattering cross section of the hydrogen atoms whose dynamics lying in the nanosecond time-and the nanometer length-scales is accessible to the spectrometer (See Methods in SM). In particular, the predominant contribution is due to the selfdiffusive dynamics of an average protein in the bacterial cytoplasm, as more than 75% of the E. coli dry weight consists of proteins and ribosomes, themselves made up of about 50% proteins and 50% RNA by mass [START_REF] Carl Neidhardt | Physiology of the bacterial cell; a molecular approach[END_REF][START_REF] Jasnin | In vivo measurement of internal and global macromolecular motions in Escherichia coli[END_REF][START_REF] Martinez | High protein flexibility and reduced hydration water dynamics are key pressure adaptive strategies in prokaryotes[END_REF][START_REF] Divina B Anunciado | In vivo protein dynamics on the nanometer length scale and nanosecond time scale[END_REF][START_REF] Mamontov | Microscopic diffusion processes measured in living planarians[END_REF]. In addition, the number of hydrogen atoms in proteins is larger than in nucleic acids, while the remainder biomolecular components, mainly phospholipids, contribute to just about 10 % of the dry weight of the bacterium. QENS allows to directly measure the incoherent dynamic structure factor S(Q, E), which can be in turn modeled to obtain unique information on the microscopic spatial and time correlation functions [START_REF] Bée | Quasielastic neutron scattering[END_REF]. The dynamical state of the proteome was sampled at increasing temperatures starting from 276 K, where the bacteria can live and thrive, till 350 K, i.e. well above the temperature of celldeath (T CD ⇡ 323 K) [START_REF] Dill | Physical limits of cells and proteomes[END_REF]. In addition, to test the reversibility of dynamical changes, we performed a few measurements while cooling the bacteria after they underwent thermal death.

Protein dynamics at all the temperatures are well described in terms of two distinct diffusive processes arising from the dynamics of the average protein in the E. coli cytoplasm (See Fig. 4.19A), whose narrow and broad Half Width-Half Maximum (HWHM) of the signal γ G and γ L correspond, respectively, to the long and the fast characteristic timescales of the probed global (G) and local (L) motions. The wavevector dependence of both γ G and γ L follows the Singwi and Sjölander jump-diffusion model [START_REF] Singwi | Diffusive motions in water and cold neutron scattering[END_REF], as shown in Fig. 4.19B and 4.19C, where the scatterers alternate oscillatory motions around their equilibrium positions for a certain residence time τ and diffusive motions between two equilibrium positions (for further details see section 4.2.2). This is consistent with results of previous in vivo QENS experiments on bacteria [START_REF] Jasnin | In vivo measurement of internal and global macromolecular motions in Escherichia coli[END_REF][START_REF] Martinez | High protein flexibility and reduced hydration water dynamics are key pressure adaptive strategies in prokaryotes[END_REF][START_REF] Mamontov | Microscopic diffusion processes measured in living planarians[END_REF]. We interpret the γ G component as due to global diffusive dynamics of an average protein while caged in the crowded E. coli cytoplasm. The faster γ L contribution arises from the local internal and inter-domain dynamics occurring within proteins [START_REF] Roosen-Runge | Protein selfdiffusion in crowded solutions[END_REF]. The apparent diffusion coefficient D G , which combines both the translational and rotational motions of proteins [START_REF] Pérez | Evolution of the internal dynamics of two globular proteins from dry powder to solution[END_REF] (see section 4.2.2) exhibits a dramatic non-reversible reduction in proximity of the cell-death temperature (Fig. 4.19D). In protein crowded solutions a similar slow-down of the diffusive dynamics was ascribed to the gelation of the system induced by the protein unfolding [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF][START_REF] Grimaldo | Hierarchical molecular dynamics of bovine serum albumin in concentrated aqueous solution below and above thermal denaturation[END_REF][START_REF] Matsarskaia | Evolution of the structure and dynamics of bovine serum albumin induced by thermal denaturation[END_REF]. The transition to a gel-like phase is also supported by the slight increase of τ G above T CD (inset of Fig. 1D), which suggests an increasing localization of proteins in cage-like structures.

Further, the local dynamics is very sensitive to temperature change, but in this case the diffusion coefficient D L shows a significant increase above T CD , while τ L is nearly constant at ⇡ 180ps (Fig. 1E and its inset). This suggests that the changes in the local dynamics are due to a variation in the extent of the explored spatial region. This is confirmed by the increase of the distance between two jumps r L , that displays a sudden break toward higher values after T CD (Fig. 4.19G). Interestingly, also the number of atoms too slow to be visible by QENS, p L , shows a similar increasing trend above T CD , as seen in Fig. 4.19F, when temperature is rising. However, comparing the fraction of H-atoms that are fixed with those that diffuse in a sphere or undergo a 3-jump diffusion (see Fig. 4.20), we find that the percentage of H-atoms that are diffusing in a spherical volume is constant in temperature. Therefore, the changes in the fraction of fixed atoms are due to the variation of the fraction of H-atoms undergoing the 3-jump diffusion, only. This behavior is consistent with the progressive unfolding of a part of the E. coli proteome. In fact, moving subgroups of denatured proteins can access larger spatial regions, and at the same time their number gradually decreases because of the growing interactions of the gel-like system formed in the cytoplasm. The rate at which protein groups locally explore their environment seems to be reversible across the thermal death, as we can see from the trend of both D L and τ L in Fig. 4.19E, confirming previous results on crowded protein solutions [START_REF] Grimaldo | Hierarchical molecular dynamics of bovine serum albumin in concentrated aqueous solution below and above thermal denaturation[END_REF]. On the other hand, the characteristic lengths of the local dynamics, i.e. r L and p L , show an irreversible trend.

Simulations show how protein unfolding slows down protein diffusion. To examine how protein motions are affected by heating and thermal denaturation of parts of the E. coli interior, we performed multi-scale molecular dynamics simulations. Previous works studied protein diffusion in crowded environments using coarse-grained (CG) [START_REF] Ando | Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion[END_REF][START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF] and all-atom [START_REF] Yu | Biomolecular interactions modulate macromolecular structure and dynamics in atomistic model of a bacterial cytoplasm[END_REF][START_REF] Sören Von Bülow | Dynamic cluster formation determines viscosity and diffusion in dense protein solutions[END_REF] MD simulations. However, while simplified CG models help sample long time-scales, they either lack solvent mediated correlations [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF] -essential for describing the protein mobility -or rely on oversimplified description of the protein lacking chemical features [START_REF] Ando | Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion[END_REF]. Here, we combined two levels of description: a CG model for protein (OPEP) owning residue level chemical resolution [START_REF] Sterpone | The OPEP protein model: from single molecules, amyloid formation, crowding and hydrodynamics to DNA/RNA systems[END_REF] to sample the local structure of the crowded solution and an all-atom description to subsequently explore the diffusion of proteins in sub-volumes obtained from the coarse-grained simulation. In our approach we used the lattice Boltzmann molecular dynamics (LBMD) technique that allows to include naturally hydrodynamic interactions in the simulations of implicit solvent molecular models, and already successfully applied to complex biological situations like protein crowded solutions [START_REF] Timr | The Unfolding Journey of Superoxide Dismutase 1 Barrels under Crowding: Atomistic Simulations Shed Light on Intermediate States and Their Interactions with Crowders[END_REF] and multi-scale protein aggregation [START_REF] Chiriccotto | Multiscale Aggregation of the Amyloid AB16-22 Peptide: From Disordered Coagulation and Lateral Branching to Amorphous Prefibrils[END_REF]. The 4.3 µs coarse-grained LBMD [START_REF] Sterpone | Protein simulations in fluids: Coupling the OPEP coarse-grained force field with hydrodynamics[END_REF] simulation contained 197 proteins of 35 different species, mimicking the protein composition of the E. coli cytoplasm [START_REF] Sean | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF] (see Fig. 4.21C for a pictorial representation of the system, and section 4.2.4 for more details on the protein composition). The length of the trajectory allowed significant reshuffling of the initial positions of the proteins and thus also exploration of different geometries of the crowded system; in fact, only 35% of the initial interaction partners were also in contact at the end of the trajectory (Fig. 4.22). The reshuffling is kinetically meaningful since the computed diffusivity for each molecular species is very close to its experimental estimate. In Fig. 4.21A we report the distribution of the translational diffusion coefficients computed for each protein. The coefficients vary between 0.5 and 7

Å2 /ns depending on the molecular weight (Table 4.4), values in excellent agreement with what reported experimentally [START_REF] Kalwarczyk | Biologistics-Diffusion Coefficients for Complete Proteome of Escherichia coli[END_REF], see blue solid lines in the graphs. From different frames of this large-scale trajectory, we selected five sub-volumes of the whole system containing 11-20 proteins so to reflect the protein composition and concentration heterogeneity of the cytoplasm (Figure 4.21B). Each sub-box was converted to the all-atom resolution (Figure 4.21C) and exposed to a sequence of production simulations (⇠ 100 ns per run) at increasing temperatures to investigate protein translational diffusion coefficients, probed in the 0.3-5 ns regime. Subsequently, for these atomistic systems we repeated this heating protocol with the same box after exposing it to a simulated rapid heat shock (see section 4.2.7), serving to completely unfold all the proteins within the time scale accessible to the simulation.

Our simulations showed a strong decrease in the average translational diffusion coefficient upon unfolding (Figure 4.21D). The observed decrease quantitatively agreed for two distinct force fields, Amber99SB-disp [START_REF] Robustelli | Developing a molecular dynamics force field for both folded and disordered protein states[END_REF] and CHARMM36m [START_REF] Huang | CHARMM36m: an improved force field for folded and intrinsically disordered proteins[END_REF], which we used to model the proteins (see Figure 4.25). In addition, for both folded and unfolded proteins the diffusion coefficient scales linearly with temperature. A similar slow-down upon unfolding has been observed in previous theoretical and experimental studies of monocomponent protein solutions [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF][START_REF] Grimaldo | Hierarchical molecular dynamics of bovine serum albumin in concentrated aqueous solution below and above thermal denaturation[END_REF][START_REF] Matsarskaia | Evolution of the structure and dynamics of bovine serum albumin induced by thermal denaturation[END_REF].

Unfolding changes protein interactions. A detailed look at the protein-protein interactions reveals the cause of the diffusion slow-down. We consider as an example one of the simulated atomistic systems of protein concentration 288 g/L. While each protein had, on average, 2.1 interaction partners in the folded system, this number increased to 8.6 when the proteins are all unfolded (see Figure 4.21E). As a consequence, the average number of atom-atom contacts per protein rose by a factor of 25 upon unfolding. Among the different contact types, this increase was the strongest for "hydrophobic" contacts (i.e., between non-polar atoms), which were enhanced by a factor of 41 (see Figure 4.21F). Recently, the diffusion slow-down in crowded solutions of globular proteins at intermediate concentrations ( 200 g/L) has been linked to the formation of transient protein clusters [START_REF] Nawrocki | Slow-Down in Diffusion in Crowded Protein Solutions Correlates with Transient Cluster Formation[END_REF][START_REF] Sören Von Bülow | Dynamic cluster formation determines viscosity and diffusion in dense protein solutions[END_REF]. Indeed, we observed that in folded systems with protein concentrations below 200 g/L, the proteins were organized in 1-3 clusters, with a few remaining protein molecules floating freely in solution. At higher concentrations, the folded proteins formed a single large cluster most of the time, while maintaining a degree of dynamical exchange with the bulk (see Figure 4.23). On the contrary, regardless of the protein concentration, unfolding led to the formation of a single cluster, encompassing all the proteins and creating a stable network. The enhanced stickiness of the unfolded proteins was reflected by a strong rise in the viscosity of the crowded protein solutions. For a 288 g/L protein system at T = 300 K, the viscosity increased upon unfolding from η = 8 ± 1 mPa•s to 60 ± 8 mPa•s with a99SB-disp and from η = 11 ± 3 mPa•s to as much as 360 ± 70 mPa•s with CHARMM36m (Fig. 4.15). The substantial relative increase in viscosity caused by unfolding persisted even at the highest temperatures, (with η going up from 3.4 ± 0.8 to 18 ± 4 mPa•s for Amber99SB-disp and from 3.7 ± 0.3 to 72 ± 19 mPa•s for CHARMM36m at T = 380 K). Our results are consistent with the trend observed in previous experimental studies [START_REF] Ml | The effect of denaturation on the viscosity of protein systems[END_REF][START_REF] Choi | Microfluidic rheometer for characterization of protein unfolding and aggregation in microflows[END_REF][START_REF] Muñoz | Detecting protein folding by thermal fluctuations of microcantilevers[END_REF], reporting an increase in the viscosity of protein solutions undergoing thermal denaturation. Owing to the significantly lower protein concentrations (below 100 g/L) that were examined in those measurements, the unfolded protein viscosities (⇠1-40 mPa•s) are lower than our computational estimates. On the other hand, the viscosities of our unfolded systems are smaller than those reported for concentrated antibody solutions [START_REF] William | Viscosity behavior of high-concentration protein mixtures[END_REF], reaching up to 1700 mPa•s, and several phase-separated biomolecular condensates [START_REF] Ghosh | Shear relaxation governs fusion dynamics of biomolecular condensates[END_REF]. In fact, given the very slow decay of the pressure autocorrelation function (Fig. 4.13), the viscosity values obtained for the unfolded systems likely represent a lower bound for the actual viscosity. Simulations show concentration dependence of diffusion slow-down. In line with previous results for crowded protein solutions [START_REF] Roosen-Runge | Protein selfdiffusion in crowded solutions[END_REF][START_REF] Grimaldo | Diffusion and dynamics of γ-globulin in crowded aqueous solutions[END_REF] and cell lysates [START_REF] Grimaldo | Protein short-time diffusion in a naturally crowded environment[END_REF][START_REF] Golub | The effect of crowding on protein stability, rigidity, and high pressure sensitivity in whole cells[END_REF], the translational diffusion coefficients extracted from the atomistic simulations decreased with increasing protein concentration (see Figure 4.24A), and the effect is stronger for the unfolded systems. In both folded and unfolded cases, the sharpest drop occurred below 200 g/L, whereas the concentration sensitivity was predicted to be weaker in the concentration range of 200-400 g/L, characteristic of the E. coli cytoplasm [START_REF] John | Macromolecular crowding: an important but neglected aspect of the intracellular environment[END_REF]. This indicates that variations in the macromolecular concentrations inside the E. coli proteome should not affected strongly the comparison of simulations with experiments. Even partial proteome unfolding causes strong diffusion slow-down. The simulations clearly prove that massive unfolding leads to a strong diffusion slow-down. However, only a fraction of the proteome might be unfolded near the cell-death temperature as suggested by recent experiments [START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF][START_REF] Mateus | Thermal Proteome Profiling in Bacteria: Probing Protein State in vivo[END_REF]. To explore the effect of such partial unfolding on the overall protein diffusivity, we performed additional temperature scans for a selected atomistic system (288 g/L) with a varying fraction r u of unfolded proteins (25%, 50%, and 75%; see Figure 4.24C). In each case we completely unfolded the chosen amount of proteins while leaving the remainder fully folded. Since the Dill model predicts similar denaturation temperatures (⇡ 328.5K ± 1.5K) for all the proteins contained in the system, we selected the proteins to be unfolded randomly, trying to maximize the species' heterogeneity of the chosen sub-set. Our simulations revealed that with the increasing unfolded content, the overall translational diffusion coefficients quickly approached those calculated for the fully unfolded system (see Figure 4.24D where we report data for Amber99SB-disp force field), this effect is even stronger when using the CHARMM36m force field (Figure 4.25). This finding demonstrates that the translational diffusion coefficient, experiencing a sharp drop already for small values of r u , is a non-linear function of r u . To explain the reason why D t shows such a rapid drop, we analysed separately the diffusion coefficients of folded and unfolded proteins in the intermediate boxes. We found that the D t of folded proteins decreased by more than 50% already for the smallest fraction of unfolded proteins (see Figure 4.24E). Thus, the presence of even a small amount of unfolded proteins is able to strongly affect the diffusion of the remaining folded proteins. On the other hand, the diffusion coefficient of the unfolded proteins shows a more gradual decrease (Figure 4.24F). On a quantitative point of view, the translational diffusion coefficient in the partially unfolded proteome can be expressed as

D t =(1 a u ) D ( f ) t + a u D (u)
t , where a u is an "apparent" unfolded fraction, weighting the diffusion coefficient of the fully folded D ( f ) t and the fully unfolded D (u) t systems. The non-linear dependence of a u of the actual proteome unfolded fraction r u describes the diffusion slow-down due to the unfolding. We found that this dependence can be fitted with the power law a u = r p u (see the insets in Figures 4.24D and 4.25C), with the exponent p equal to 0.411 ± 0.026 for a99SB-disp and 0.142 ± 0.012 for CHARMM36m. Unlike the translational diffusion coefficient, the dependence of viscosity on the unfolded content appeared to be nearly linear, within the limits posed by a slow convergence (Figure 4.26). As we describe in section 4.4, a linear unfolding dependence of viscosity combined with the simplified Stokes-Einstein model results in a non-linear dependence of a u on r u that is similar to the one recovered directly from simulations (Figure 4.27).

Estimation of the folding state of the cytoplasm and validation of the results. By combining the QENS and MD simulations results, we estimated the temperature dependence of the fraction of unfolded proteins in the E. coli cytoplasm.

For this purpose, and inspired by previous work on protein crowded solutions [START_REF] Grimaldo | Hierarchical molecular dynamics of bovine serum albumin in concentrated aqueous solution below and above thermal denaturation[END_REF], we first consider the temperature dependent profile of the apparent diffusion coefficient measured from QENS experiments (Fig. 4.19D), and we perform a fit using an empirical relation that combines the diffusivity of proteins in the folded and unfolded states,

D(T )=D ( f ) • [1 a QENS u ]+D (u) • a QENS u (T )
, where the function a QENS u , that measures the changes in the global diffusion, is equivalent to the apparent fraction of unfolded proteins in the system as defined in the previous paragraph. In the second step we exploited the empirical relationship derived from simulations between the real fraction of unfolded proteins in the solution, r u , with its apparent counterpart, a u . We applied at each temperature the empirical relation to the experimentally derived apparent fraction of unfolded proteins 1

, r QENS u (T )= h a QENS u (T ) i 1/p
. With this tool in hands, we obtained the experimentally derived fraction of unfolded proteins as a function of the temperature, r QENS u (T ), see Fig. 4.28A. This quantity is force field dependent via the exponent parameter p, and can be now compared with the equivalent function theoretically derived by Dill and coworkers [START_REF] Dill | Physical limits of cells and proteomes[END_REF], see Fig. 4.28A.

Strikingly, at variance with the proteome catastrophe scenario derived theoretically [START_REF] Dill | Physical limits of cells and proteomes[END_REF], the experimentally derived r QENS u s show a very slow increase of the fraction of unfolded proteins as temperature crosses T CD = 323.15K . For sake of clarity, a few degrees above the cell death temperature, the Dill's model predicts more than 50% of the proteins unfolded, while the experimentally derived unfolded fraction is, on average, less than 15% (see Table 4.10). This result is in line with recent experiments on proteome thermal stability [START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF][START_REF] Mateus | Thermal Proteome Profiling in Bacteria: Probing Protein State in vivo[END_REF]. It is worth to recall that our approach is based on two assumptions (see section 4.2). First, we hypothesized that the main contribution to the QENS signal comes from an average protein. This is a quite reasonable approximation, as ribosomes, DNAs and RNAs are too massive and consequently slow to be detected, while the remainder biomolecular components, mainly phospholipids, contribute to just about the 10% of the dry weight of the bacterium [START_REF] Carl Neidhardt | Physiology of the bacterial cell; a molecular approach[END_REF]. In addition, in our simulations we represented a very simplified version of the E. coli cytoplasm, composed of just a small subset of proteins. This simplified representation, however, is able to catch the main dynamic features of the system we investigated, thus strengthening the picture we propose. (see the Discussion for more details).

An additional support to our findings comes from the fact that the experimental apparent diffusion coefficient D G of the average protein can be correctly described from the simulations just starting from the estimates of the apparent fraction of unfolded proteins and of D ( f ,u) G . The curves are reported in Fig. 4.28B, showing a very good agreement with the experimental data, thus endorsing the assumptions we made.

Reproducing the growth rate of E. coli. The stability curve of the proteome is used now to reconstruct the growth-rate curve of E. coli. We follow in spirit the approach described by Dill in [START_REF] Dill | Physical limits of cells and proteomes[END_REF] where the growth-rate (g(T )) is related to the temperature dependent fraction of unfolded protein via an Arrhenius reaction rate term, g(T )=g 0 e ∆H † ∏ Γ i=1 r i (T ); where g 0 is an intrinsic growth-rate parameter, ∆H † is the dominant activation barrier, Γ is the number of essential proteins for the bacterium growth, and r i is the temperature dependent fraction of unfolding for the protein species i. In our approach, r i is replaced by the average estimated r QENS u . When using r QENS u from a99SB-disp we obtained an excellent fit, and values for ∆H † ' 45 kJ/mol and

Connection between the apparent and the real unfolded fractions

As we detailed in the subsection "Neutron scattering data analysis", we found that D G measured by QENS experiments can be described as:

D G (T )=(T • a 1 + b 1 ) • [1 a QENS u (T )] + (T • a 2 + b 2 ) • a QENS u (T ) (4.27)
enough to ensure that the folding state of the proteins was not affected. This allowed us to verify that if the folding state remains constant, the diffusion coefficient of the average protein scales linearly with the temperature (Figure 3D and Comparing the diffusion coefficients calculated for a system simulated in different folding states (Figure 3D and4.25C), we found a relationship between the amount of unfolded proteins and the consequent diffusion slowdown-eq. (4.30). Finally, with the inverse relation eq. (4.31) and the QENS measurement of the slowdown, we were able to estimate the real fraction of unfolded proteins at different temperatures in the E. coli cytoplasm through the study of its dynamical state (see Figure 4A).

Origins of the nonlinearity of the transfer function

The observation that the translation diffusion coefficient D t cannot be expressed simply as

D t (r u )=(1 r u )D ( f ) t + r u D (u) t (4.32)
where r u is the unfolded fraction and

D ( f ) t and D (u)
t are the translational diffusion coefficients corresponding to a fully folded-and a fully unfolded system, respectively, can be rationalized in the following way. If we assume the validity of the Stokes-Einstein relation, the diffusion coefficients

D ( f ) t and D (u)
t , forming two limiting cases, can be written as

D ( f ) t = k B T 6πη ( f ) R ( f ) (4.33)
and

D (u) t = k B T 6πη (u) R (u) (4.34)
Here, R ( f ) and R (u) are protein hydrodynamic radii in the folded and unfolded states, respectively, and η ( f ) and η (u) denote the viscosities of the fully folded and fully un- folded volumes.

However, in intermediate boxes with an unfolded fraction r u , η(r u ) is neither identical to η ( f ) nor to η (u) . Consequently, for the diffusion of folded/unfolded proteins in the intermediate boxes, we can expect

D ( f ) t (r u )= k B T 6πη(r u )R ( f ) (4.35)
and

D (u) t (r u )= k B T 6πη(r u )R (u) (4.36)
where the viscosity is a function of r u . Therefore, rather than assuming Eq. 4.32, one should take into account the effect of the intermediate viscosity:

D t (r u )=(1 r u )D ( f ) t (r u )+r u D (u) t (r u )=(1 r u ) η ( f ) η(r u ) D ( f ) t + r u η (u) η(r u ) D (u) t (4.37)
If we consider a linear dependence of η on r u ,

η(r u )=η ( f ) + r u (η (u) η ( f ) ) (4.38)
with values of η ( f ) and η (u) estimated using simulations, we obtain a similar depen- dence of a u on r u as the one based on the actual diffusion coefficients determined from simulations (see Figure 4.27).

Further Results

scientific results obtained during the Ph.D. studies concerning the effects of protein-ligand complexation on protein dynamics In this chapter, we present a study comparing atomic motional amplitudes in calcium rich and depleted alpha-lactalbumin. The investigations were performed by elastic incoherent neutron scattering (EINS) and molecular dynamics (MD) simulations. As the variations were expected to be very small, three different hydration levels and timescales (instrumental resolutions) were measured. In addition, we used two models to extract the mean square displacements (MSDs) from the EINS data, one taking into account the motional heterogeneity of the MSD. At a timescale of several nanoseconds, small differences in the amplitudes between the calcium enriched and depleted alpha-lactalbumin are visible, whereas at lower timescales no changes can be concluded within the statistics. The results are compared to MD simulations at 280 and 300 K by extracting the MSDs of the trajectories in two separate ways: first by direct calculation, and second by a virtual neutron experiment using the same models as for the experimental data. We show that the simulated data give qualitatively similar results as the experimental data but quantitatively there are differences. Furthermore, the distribution of the MSDs in the simulations suggests that the inclusion of heterogeneity is reasonable for alpha-lactalbumin, but a bi-or trimodal approach may be sufficient.

Introduction

A number of proteins have the ability of binding ions that may lead to changes in the protein's structure and dynamics at the atomic scale, and subsequently, may affect their functionality. Recent studies have shown that for example, in the case of enzymes the presence of small inhibitors might influence the dynamics in a measurable way compared to the dynamics of the wild type form [START_REF] Nitta | The evolution of lysozyme and αlactalbumin[END_REF][START_REF] Trapp | Correlation of the dynamics of native human acetylcholinesterase and its inhibited huperzine A counterpart from sub-picoseconds to nanoseconds[END_REF][START_REF] Peters | Dynamics of human acetylcholinesterase bound to non-covalent and covalent inhibitors shedding light on changes to the water network structure[END_REF][START_REF] Saouessi | Asymptotic analysis of quasielastic neutron scattering data from human acetylcholinesterase reveals subtle dynamical changes upon ligand binding[END_REF][START_REF] Saouessi | Frequency domain modeling of quasielastic neutron scattering from hydrated protein powders: Application to free and inhibited human acetylcholinesterase[END_REF] and single point genetic mutation in proteins can affect collective density fluctuations in hydrating water [START_REF] Russo | Water collective dynamics in whole photosynthetic green algae as affected by protein single mutation[END_REF]. Another case is alpha-lactalbumin (α-La), the major whey protein found in the milk of all mammals. It is a simple Ca 2+ binding milk protein and has a significant role in biosynthesis of lactose in the lactating mammary gland. Together with the enzyme β -1,4-galactosyltransferase (β 4GalT) it forms a complex which is responsible for the lactose synthase, i.e., transforming galactose and glucose into lactose. It strongly binds the cation Ca 2+ and results in changes in the tertiary structure of the protein (see Fig. 5.1). Besides Ca 2+ , the binding site can also bind Mg 2+ , Mn 2+ , Na + , or K + , which induce similar but smaller structural changes than Ca 2+ . However, the corresponding binding constants are much lower except in the case of Mn 2+ . In general, the binding of a cation stabilizes α-La and increases its thermal denaturation temperature. Furthermore, recently Shinozaki and Iwaoka [START_REF] Shinozaki | Effects of metal ions, temperature, and a denaturant on the oxidative folding pathways of bovine α-lactalbumin[END_REF] showed that Ca 2+ and Mn 2+ accelerates folding to the native form of α-La, an effect not seen with the other cations. α-La can also bind Zn 2+ at several other distinct binding sites, but results in a decrease in the stability of α-La bound to Ca 2+ . The apo form refers to α-La which is not bound to Ca 2+ . Owing to the characteristics described above, α-La is often used as a simple model for Ca 2+ binding proteins.

In addition to structural changes, the binding of α-La to Ca 2+ may also generate structural rearrangements capable of influencing locally molecular dynamics and therefore varying the functionality of the protein. The task of probing such small effects is not easy and a sophisticated approach is required. Incoherent neutron scattering is a technique used to probe atomic and molecular dynamics on timescales of pico-to nanoseconds, and when combined with molecular dynamics (MD) simulations forms a powerful partnership, and could indeed offer unique insights into changes occurring in the atomic scale (amplitudes of about a few angstroms). Despite the two techniques accessing very similar times and dimensions, they are not always in full quantitative agreement [START_REF] Bachir Aoun | Direct comparison of elastic incoherent neutron scattering experiments with molecular dynamics simulations of DMPC phase transitions[END_REF]. To better understand the reasons behind such disagreements, we choose to make a detailed study of the dynamics of α-La, which is commercially available, to do exhaustive neutron experiments, and sufficiently small to permit accurate simulations. We study the dynamics in both the Ca rich and depleted forms, which from hereafter will be referred to as α-La ca and α-La dep , respectively.

The intensities measured using incoherent neutron scattering experiments are commonly used to extract mean square displacements (MSDs) of the protons within the protein, almost exclusively by assuming a harmonic approximation of all possible dynamical contributions, i.e., the Gaussian approximation [START_REF] Azms Rahman | Theory of slow neutron scattering by liquids[END_REF][START_REF] Squires | Introduction to the theory of thermal neutron scattering[END_REF]. Furthermore, it is common to combine results from different neutron spectrometers, since they cover different timescales and length scales. Recently, we applied models that go beyond this approximation and include dynamic heterogeneity, to be able to fully exploit a wider instrumental spatial window [START_REF] Zeller | Analysis of elastic incoherent neutron scattering data beyond the Gaussian approximation[END_REF]. In this work, we combine data from three neutron spectrometers which access different timescales (have different energy resolutions) and length scales (have different momentum transfer coverage). We apply a few models, including the commonly used Gaussian approximation, to the data to investigate to what extent they help to disentangle small effects on the dynamics and make a comparison to results from MD simulations.

Experimental Section

Sample preparation for neutron experiments

All experiments described use bovine alpha-lactalbumin (α-La), either in its natural form with Ca 2+ (α-La ca ) or Ca 2+ depleted α-La dep . The protein was purchased from Sigma-Aldrich in lyophilized powder form. Three different hydration levels were prepared for each batch and protein type, hydrated with heavy water, D 2 O. This is so the neutron signal is dominated by the incoherent scattering from the protons in the protein (owing to the large incoherent neutron cross section of hydrogen compared to deuterium or other atoms constituting the protein structure [START_REF] Varley | Neutron scattering lengths and cross sections[END_REF]). The hydration level was determined from the difference in mass with and without D 2 O and is defined as h = grams D 2 O/grams dry protein. The different levels of hydration were h ⇡ 0 (dry), h ⇡ 0.4, and h ⇡ 0.8. The dry lyophilized sample represents the case where only harmonic motions are present up to room temperature and 0.4h corresponds to around a hydration level of one or two layers of water on the protein surface [START_REF] Perticaroli | Description of hydration water in protein (green fluorescent protein) solution[END_REF], which is sufficient to allow for localized dynamical motions. Finally, 0.8h represents a gel state close to full hydration. The purchased lyophilized protein powder was dried for at least 24 h, after which it was weighed and then loaded into flat aluminum sample holders (standard for neutron spectroscopy experiments) and vacuum sealed with indium wire. The 0.4h and 0.8h samples were also dried before hydration, and then left in a D 2 O rich environment to uptake the water. The samples were weighed periodically until they achieved the desired uptake of D 2 O and then sealed with indium in similar flat aluminum holders. Masses of around 100 mg were used to obtain around 10% scattering and ensure there was no significant multiple scattering.

Simulation setup

MD simulations of hydrated protein powder were used representing the interactions of proteins with a small amount of water, to be able to compare them directly with experimental data. The simulations were started from two different protein structures which can be found in the Protein Data Bank (PDB) [START_REF] Helen M Berman | The protein data bank[END_REF]: (1) bovine α-La with calcium (α-La ca ), PDB ID: 1F6S and (2) bovine α-La without calcium (α-La dep ), PDB ID: 1F6R. Both structures were published by Chrysina et al. [START_REF] Evangelia D Chrysina | Crystal structures of apo-and holo-bovine α-lactalbumin at 2.2-Åresolution reveal an effect of calcium on inter-lobe interactions[END_REF]. Each PDB structure consists of six distinct α-La proteins (= chains), allowing to calculate an average dynamics of a single α-La chain. As a matter of fact, Tarek and Tobias [START_REF] Tarek | The dynamics of protein hydration water: a quantitative comparison of molecular dynamics simulations and neutronscattering experiments[END_REF] pointed out that a single protein covered by a shell of water is not sufficient to describe a powder protein by simulations. Instead, a crystal composed by two proteins or more resulted in a realistic model to reproduce neutron scattering data. This is the reason why, in the present case, we used six chains of proteins placed in each box (see Fig. 5.

2).

The protein molecule was centered in a cubic box of size 8.39 nm at first, with the CHARMM27 force field [START_REF] Jr | All-atom empirical potential for molecular modeling and dynamics studies of proteins[END_REF][START_REF] Jr | Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations[END_REF], and the TIP4P-EW water model [START_REF] Hans W Horn | Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew[END_REF], using GROMACS 5.0.7 (GPU version) as the MD engine [START_REF] Van | GROMACS: fast, flexible, and free[END_REF][START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF]. The boxes were filled with water molecules to start with, which were then deleted (starting from the outside) until the number of water molecules around the protein met the desired hydration level h. The box for a hydration level of 0.4h contained 1824 (dep) / 1834 (Ca) water molecules and 232 (both) for the dry system (0.05h). All systems were electrically neutralized by adding NaCl. Van der Waals interaction was truncated at 1.2 nm with the Lennard-Jones potential switched to zero gradually at 1.0 nm. A particle mesh Ewald [START_REF] Essmann | A smooth particle mesh Ewald method[END_REF] with a Coulomb cutoff of 1.2 nm was used to calculate electrostatic interaction. All bonds involving hydrogen atoms were constrained with the LINCS [START_REF] Hess | LINCS: a linear constraint solver for molecular simulations[END_REF] algorithm. The systems were firstly energetically minimized using steepest descent steps with a maximum force of 10.0 kJ mol 1 nm 1 and a maximum of 50000 steps. Then they were equilibrated in the NVT ensemble at T = 280 K (and 300K) for 300 ps and in the NPT ensemble at p = 1 bar for 50 ns, with a 0.5-fs time step to slowly release the unreasonable atom contact and suppress vacuum. The temperature coupling was performed using the velocity-rescale algorithm with a coupling time of 1 ps [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF]. The pressure coupling was performed using the Parrinello-Rahman algorithm with a coupling time of 1 ps [START_REF] Melchionna | Hoover NPT dynamics for systems varying in shape and size[END_REF]. The production MD simulations for hydration level 0.4h were conducted in the NPT ensemble for 100 ns, with a 2-fs time step, while those of the dry systems were conducted for 500 ns. Only the last 20-ns trajectories recorded at every 2 ps were used for the analysis. For such a dense system, the global translation and rotation of the protein molecules was strongly suppressed [START_REF] Jh Roh | Onsets of anharmonicity in protein dynamics[END_REF][START_REF] Tan | Gradual crossover from subdiffusion to normal diffusion: a many-body effect in protein surface water[END_REF]. A visualization of the difference in box size and hydration level is shown in Fig. 5.2.

Neutron Scattering

Data were collected on three neutron spectrometers, all so-called inverted geometry spectrometers, covering a wide temporal range, namely, OSIRIS [START_REF] Mark | Spectroscopic characteristics of the OSIRIS near-backscattering crystal analyser spectrometer on the ISIS pulsed neutron source[END_REF] at the ISIS Neutron and Muon Facility, UK; IN13 [START_REF] Francesca | IN13 backscattering spectrometer at ILL: looking for motions in biological macromolecules and organisms[END_REF] at the Institut Laue Langevin, (ILL), Grenoble, France; and SPHERES [START_REF] Wuttke | Jülich's high-flux neutron backscattering spectrometer at FRM II[END_REF] at the MLZ Munich reactor in Germany. The data collected at the ILL can be found under the DOIs in Refs. [START_REF] Zeller | Complete data set of different hydration levels of bovine alpha-Lactalbumin[END_REF] and [START_REF] Zeller | Testing the Validity of Current Models to Describe the Protein Dynamics from EFWS data[END_REF]. OSIRIS and SPHERES use crystal analyzers that reflect cold neutrons (λ of 6.27 and 6.66 Å, respectively) allowing access to a momentum transfer range, q-range, up to 1.8 Å 1 , whereas IN13 uses a thermal neutron crystal analyzer (λ of 2.23 Å) which opens up the accessible q-range to 4.9 Å 1 . This permits us to probe dynamics occurring in a variety of length scales, where distance d = 2π/q. In addition, the three instruments differ in energy resolutions allowing access to motions from a few picoseconds to a few nanoseconds. Specifically, they are 25, 8, and 0.7 µeV, for OSIRIS, IN13, and SPHERES, respectively. Transmission values for all samples were measured on IN13 to be above 90% so that multiple scattering effects were not taken into consideration for the data treatment. The initial data reduction was done with LAMP [START_REF] Richard | Analysis and Visualisation of Neutron-Scattering Data[END_REF] for IN13, SLAW for SPHERES, and MANTID [START_REF] Arnold | Mantid-Data analysis and visualization package for neutron scattering and µ SR experiments[END_REF] for OSIRIS. Slab can corrections for a flat sample holder and normalizations providing the relative detector efficiency and the instrumental resolution were done with LAMP for the samples measured on IN13 and SPHERES. The measurements on OSIRIS were corrected using the empty sample holder and normalized in MANTID. All intensity normalizations were done with the lowest available temperature data of each scan. Therefore, the difference between the slab correction algorithm and the subtraction of the empty sample holder alone are negligible.

Incoherent neutron scattering measurements give access to the elastic incoherent structure factor (EISF), S, which is a function of the momentum transfer q at the elastic line, where the energy transfer hω that occurs between the neutrons and the scattering atoms (mostly hydrogen), as a result of the scattering event, is approximately zero. The most commonly used approach to analyze this intensity is to assume that the atomic nuclei undergo harmonic motions around their equilibrium positions [START_REF] Azms Rahman | Theory of slow neutron scattering by liquids[END_REF] and thus fit the data to the so-called Gaussian approximation (GA). The intensity can then be expressed as

S(q, 0 ± ∆E; ⌦ r 2 ↵ ) ⇡ S 0 exp q 2 ⌦ r 2 ↵ GA 3 ! (5.1)
where ∆E corresponds to the instrumental energy resolution. From this expression, values for the static mean square displacements of the atoms, ⌦ r 2 ↵ GA , are obtained at each temperature point measured, by fitting the slope of the logarithm of the scattered intensities plotted vs q 2 according to

⌦ r 2 ↵ GA ⇡ 3 ∂ ln S(q, 0 ± ∆E; ⌦ r 2 ↵ ) ∂ q 2 (5.2)
The Gaussian approximation is strictly valid for q ! 0, and it holds up to q 2 ⌦ r 2 ↵ GA ⇡ 1, restricting the q-range that can be used for this type of analysis considerably.

A model that imposes no constraints on the q range is that developed by Kneller and Hinsen [START_REF] Gerald | Quantitative model for the heterogeneity of atomic position fluctuations in proteins: a simulation study[END_REF] and applied to experimental data by Peters and Kneller [START_REF] Peters | Motional heterogeneity in human acetylcholinesterase revealed by a non-Gaussian model for elastic incoherent neutron scattering[END_REF]. It differs from the GA in that it takes into account motional heterogeneity of the amino acid side chains and their environment, compared to the Gaussian approximation where only one atomic motion is representative for all hydrogens. The motional heterogeneity of the hydrogen atoms is described by a Gamma distribution and the corresponding elastic intensity can be calculated analytically as

S(q; ⌦ r 2 ↵ , β )= 1 ✓ 1 + q 2 •hr 2 i PK 3β ◆ β (5.3)
where β is a measure of the homogeneity in the atomic motions; e.g., when β ! 0 the Gaussian form is retrieved. Fits of the data give, then, access to the corresponding static mean square displacement, ⌦ r 2 ↵ PK , where PK stands for the Peters-Kneller model hereafter.

An earlier attempt to account for motional heterogeneity in modeling the EISF was suggested by Meinhold et al. [START_REF] Meinhold | Protein dynamics and stability: the distribution of atomic fluctuations in thermophilic and mesophilic dihydrofolate reductase derived using elastic incoherent neutron scattering[END_REF] by describing the mean square motional amplitudes by a Weibull distribution. However, this approach is not investigated here.

Analysis of the simulated data

Direct calculation of the MSD

The α-La ca (α-La dep ) proteins used in the simulations consist of a total of 11512 (11457) protein atoms. The number of atoms is different in the two forms, because some chains are missing some amino acids (residues) at the end of the α-La chain since they were not resolved in the PDB structure. In order to compare the simulation data to the experiment, we analyze the H atoms in the protein, which account for the majority of the scattering signal in the neutron experiments. Furthermore, in order to be consistent in the evaluation of the MSD, only H atoms which are in all chains are considered: Every single α-La consists of at least 922 H atoms which are of the same type for all α-La protein chains. Therefore, with six single α-La chains in each simulation, in total 6 ⇥ 922 = 5532 H atoms have been evaluated to calculate the MSD and thus to analyze the averaged atomic movements of the protein. The MSD of a single atom α at location r α (t) at time step t in the simulation is calculated via

MSD α (t)= D [r α (t 0 ) r α (t + t 0 )] 2 E t 0
(5.4) with h•••i t 0 being the average over all t 0 defined by the time steps of the simulation. From these individual atoms, a mean µ(t) of the MSD can be calculated. We first calculated the time average of the MSD according to Eq. (5.4) using the complete 20-ns trajectories. Further, to estimate the error of the mean of the MSD due to different conformations, the 20-ns simulations were truncated in four equally time spaced parts of 5 ns. The result of the four independent parts was then averaged to obtain a mean MSD µ i (t) and its sample standard deviation s(t) taken as an estimation of the error:

s(t)= v u u t 1 N 1 N=4 ∑ i=1 [µ i (t) µ i (t)] 2 (5.5)
Finally, we compared the MSDs obtained from a direct calculation with the ones using the fast correlation algorithm proposed by Kneller et al. [START_REF] Gerald R Kneller | nMOLDYN: a program package for a neutron scattering oriented analysis of molecular dynamics simulations[END_REF]. See Figs. 5.7 and 5.8 in the supporting information section for the results of the different checks.

In order to compare the dynamic MSD µ(t) of the simulations with the static MSD ⌦ r 2 ↵ calculated by the models, it has to be divided by 2 since the static MSD ⌦ r 2 ↵ is defined as a time independent quantity due to the confined motion resulting in [START_REF] Vural | Determination of dynamical heterogeneity from dynamic neutron scattering of proteins[END_REF][START_REF] Vural | Motional displacements in proteins: The origin of wave-vector-dependent values[END_REF]

2 ⌦ r 2 ↵ = MSD(t ! 0) (5.6)
For convenience, in the following, the MSDs obtained from the simulations will be labeled as direct MSD ∆ dir (t) with the following definition:

∆ dir (t)= 1 2 µ(t) (5.7)
The ∆ dir (t) can be directly compared to the MSDs extracted from the fits [see Eqs. (5.2) and (5.3)].

Indirect calculation of the MSD

An alternative way to compare the simulation results with the experimental data is to extract the MSD from the convolution of the instrumental resolution R(ω) with the theoretical dynamic incoherent structure factor (DISF) S inc (q, ω) calculated with the help of the simulation data. The DISF was calculated with the program MDANSE [START_REF] Goret | MDANSE: An interactive analysis environment for molecular dynamics simulations[END_REF] (v.1.1). The resolution function R(ω) for each instrument was approximated by a normalized Gaussian function with a full width at half maximum (FWHM) equivalent to the resolution of the instrument:

G(ω,t)= 1 σ res p 2π • exp ( 1 2 ✓ ω σ res ◆ 2 ) (5.8)
where FWHM = σ res p 8 ln 2 ⇡ 2.35 • σ res (5.9)

The FWHM of each instrument was obtained by matching the above-defined Gaussian function to data from vanadium which is used to measure, experimentally, the resolution of neutron spectrometers since it is an isotropic incoherent scatterer. For IN13 and OSIRIS data from a vanadium standard summed over all momentum transfers, q was used; for SPHERES, the resolution function found in the literature for the large angle detectors was used [START_REF] Wuttke | Simulation-guided optimization of small-angle analyzer geometry in the neutron backscattering spectrometer SPHERES[END_REF] (see Fig. 5.2) (Voigt profile with σ res = 0.244 µeV; γ res = 0.052µeV). The resolution functions are then convoluted with the DISF which is obtained from the simulation. For each DISF calculated with an absolute momentum transfer q m , N q = 50 q-vectors q i with a randomized direction and an absolute length of q i = q m + ∆q, with ∆q  0.05 Å 1 , are averaged. In total, the DISF is then calculated in MDANSE as

S inc (q m , ω)= 1 2π Z +∞ ∞ I inc (q m ,t) • exp(iωt) dt (5.10) I inc (q m ,t)= 1 N α N α ∑ α * 1 N q N q ∑ i exp {iq i r α (t 0 )} • exp { iq i r α (t 0 + t)} + t 0 (5.11)
where N α is the number of H atoms in the simulation and r α their location. t and t 0 are defined by the time steps of the trajectory.

From the resolution broadened DISF, S R inc (q, ∆ω), the elastic tic incoherent structure factor EISF(q m ) is computed by summing up the intensities in the range ω = ±FWHM/2 and the resulting EISF(q) is normalized by EISF(q m = 0). The obtained EISF(q) can be fitted in the same way as the experimental data to calculate the MSD. It is important to mention that for the experimental data the lowest temperature scan was used for the normalization, whereas here the value obtained at q m = 0 was taken due to the lack of a simulation at very low temperature. The models chosen to analyze the EISF(q) are the same as for the experimental data, namely, the Gaussian approximation [Eq. (5.1)] and the PK model [Eq. (5.3)], over q ranges of 0 1 Å 1 and 0 4 Å 1 , respectively.

Neutron Scattering Results

The MSDs were extracted as described above for the three instruments and according to the two models. The PK model fits data over a much wider range of q values, and thus is expected to yield additional information about the motional amplitudes. As, for instance, methyl group rotations are small motions and become particularly visible only at higher q values (> 2 Å 1 ) [START_REF] Doster | Are proteins dynamically heterogeneous? Neutron scattering analysis of hydrogen displacement distributions[END_REF][START_REF] Doster | Scaling analysis of bio-molecular dynamics derived from elastic incoherent neutron scattering experiments[END_REF][START_REF] Doster | Protein-water displacement distributions[END_REF], such treatment is expected to give a more precise description. However, it includes one more fit parameter and gives thus higher error bars for the fitting parameters. No significant differences were observed within statistical error between the dynamics (and MSD) of α-La dep and α-La ca on the timescales of OSIRIS or IN13 at any of the three hydrations. An example of this is shown for the OSIRIS data fit using the PK model, in Fig. 5.3 (top). Furthermore, both models show similar trends. Absolutely no differences are appreciable in the dry pro- teins (blue curves). At h = 0.4 (red curves) and h = 0.8 (green curves), the dynamics are almost identical except at the higher temperatures (290 -310 K), where small differences are visible and both models suggest that the α-La ca has a slightly smaller MSD, indicating less dynamics. Given the large error bars, the effect is not conclusive; however, the trend would confirm the findings of Chrysina et al. [START_REF] Evangelia D Chrysina | Crystal structures of apo-and holo-bovine α-lactalbumin at 2.2-Åresolution reveal an effect of calcium on inter-lobe interactions[END_REF], which suggest that the binding of a protein to a cation stabilizes the protein, irrespective of the hydration level.

Differences between the dynamics of the two samples are visible on the timescale of the SPHERES instrument (Figs. 3(b) and 3(c), slower dynamics up to a couple of nanoseconds). Already in the dry state, the α-La ca sample has a slightly higher MSD than the α-La dep sample above 250 K. This difference is emphasized when using the model that uses a larger q range, the PK model, suggesting that also smaller amplitudes (corresponding to higher q values) have to be included in the analysis to permit such a subtle differentiation.

At h = 0.4 no difference between the samples is observed in the PK model, and using the GA model gives a small difference where the α-La dep is more mobile than the α-La ca . This could indeed be the case, as for the highly hydrated samples (h = 0.8) the same trend is observed, and, more apparent, the MSD is larger for the α-La dep above 270 K. It is in fact the opposite behavior as for the dry, but more in line with the expected scenario of stabilization of the α-La upon binding calcium. A higher resilience of a protein upon binding of a cation indicates an increased free energy including a higher enthalpy arising from bonded interactions [START_REF] Tehei | Fast dynamics of halophilic malate dehydrogenase and BSA measured by neutron scattering under various solvent conditions influencing protein stability[END_REF]. Entropy is likely rather unchanged at the same hydration level. It appears that the sample hydrated at 0.4h presents higher dynamics between 200 and 270 K than the one hydrated at 0.8h. Similar effects were already observed for the green fluorescent protein (GFP) [START_REF] Hong | Elastic and Conformational Softness of a Globular Protein[END_REF][START_REF] Nickels | Dynamics of Protein and its Hydration Water: Neutron Scattering Studies on Fully Deuterated GFP[END_REF] and interpreted by the authors as a suppression of protein dynamics at lower temperatures by hydration water and an enhancement of it at higher temperatures. Moreover, in the 0.4h sample, the water is in a confined or glassy state so that secondary relaxations set in upon heating, whereas in the 0.8h sample where water is primarily bulk water, it is in a frozen state. The steplike increase of the MSD around 270 K for the highest hydrated sample corresponds thus to the melting of the surrounding water.

The GA model shows a more pronounced feature at the melting of ice in the higher hydrated sample compared to the PK model. The GA model covers indeed only larger length scales representing more likely the melting of the ice, whereas in the PK model, which also covers local length scales, an average of the larger and smaller length scales slightly smears out such effects.

MD Simulation Results

The MSD obtained directly and indirectly from the simulations were compared. To relate the time dependent direct MSD ∆ dir (t) to the results of the time independent indirect MSD ⌦ r 2 ↵ ind , Heisenberg's uncertainty principle was used: with τ FWHM = h/FWHM. The corresponding times for each instrument are summa- rized in Table 5.1. Figure 5.4 shows the results of the MSDs of the two different fitting models GA and PK (blue and red, respectively) and the directly calculated MSDs (black), for the two simulated sam-ples (dry and hydrated) at all three instrumental resolutions. The results obtained with the various methods to calculate the MSDs directly were so close (differences below 0.5%) that it was not possible to represent them individually in Fig. 5.4. For the dry protein, the MSD calculated via the direct method is always larger than the MSD calculated from the models, but the behavior between the simulations is the same. Assuming that the direct calculation represents a result as close as possible to the true MSD, the difference between the values of the direct method and those from the models could indicate the order of magnitude of the error introduced by using models. As anticipated, the MSD increases with increasing temperature and in fact the effect is larger on smaller timescales, i.e., at lower instrumental resolutions (IN13 and OSIRIS). The α-La ca simulations also have a slightly higher MSD for all instrumental resolutions. When comparing the different models, the GA evaluates to a higher MSD than the PK model for IN13 and OSIRIS. For SPHERES this behavior is inverted. Similar trends are observed for the hydrated protein, except for three main differences. First, the MSD of the models is much closer to the direct MSD, albeit still smaller. Secondly, the difference between the MSDs at 280 and 300 K is much larger. Thirdly, for SPHERES the MSD for α-La ca at 300 K is slightly lower than for α-La dep , which is the case for all methods considered.

The simulation also allows us to calculate the distribution of the MSDs for the protons in the protein. This is calculated following the method used by Yi et al. [START_REF] Yi | Derivation of Mean-Square Displacements for Protein Dynamics from Elastic Incoherent Neutron Scattering[END_REF] which enables an evaluation of the main contributions to the heterogeneity and of how many populations with different motions are present. Figure 5.5 shows the distribution at t = 30 ps (OSIRIS), t = 60 ps (IN13), and t = 1ns (SPHERES) for all simulations. The curve for each time was obtained by binning the individual direct MSD values in steps of 0.02 Å2 together and normalized by the total number of H atoms. The individual direct MSD values were obtained by averaging the value of the four independent slices of 5 ns for each simulation in the same way as for the direct MSD evaluation. In all simulations and for all three times t , one large peak at around 0.13 Å2 is visible (dashed vertical line). Only for the distribution at 1 ns (green) a small second peak around 1.35 Å2 is visible. The latter peak was identified in the simulations to correspond to methyl group rotations, which exist also within the IN13 data [START_REF] Peters | Motional heterogeneity in human acetylcholinesterase revealed by a non-Gaussian model for elastic incoherent neutron scattering[END_REF][START_REF] Doster | Are proteins dynamically heterogeneous? Neutron scattering analysis of hydrogen displacement distributions[END_REF], but are then retrieved at much smaller MSD values below 0.5 Å2 and cannot be separated from the motions in the main peak. For the hydrated samples the first peak is shifted slightly to higher MSD values and its peak is significantly smaller than for 30 and 60 ps. This effect is emphasized at 300 K. No significant variation is observed between the α-La dep and α-La ca samples. This seems understandable as the components forming the two samples are extremely similar.

Such analysis helps us to understand if a complete distribution of Gaussian motions is required to describe the MSD or if a bi-or trimodal approach is sufficient. According to our results, a bimodal description seems to be very reasonable, in agreement with recent works of Vural et al. [START_REF] Vural | Determination of dynamical heterogeneity from dynamic neutron scattering of proteins[END_REF] or Doster [START_REF] Doster | Are proteins dynamically heterogeneous? Neutron scattering analysis of hydrogen displacement distributions[END_REF].

Comparison of Experimental and Simulated Results

To compare the MD simulations with the experimental data, the results of the fitting models from the previous section are plotted together with the experimental results. The experimental data were collected at 5-10 K intervals, which unfortunately are not always in coincidence with the two simulation temperatures. To reduce the effects on the results, experimental MSD values were averaged over three temperature values (smoothing average) and then the MSDs between two smoothed temperature data points have been linearly interpolated. This procedure ensured that the simulated and experimental data were at the same temperature as the simulations.

As can be seen in Fig. 5.6, the MSDs extracted from the simulated data agree well with the experimental data and indicate that the different models hardly allow differentiation. Furthermore, the variations between simulated and experimental results may arise mainly from the instrumental limitations. Finally, the differences between the models are larger for the experimental data reflecting the worse statistics.

The experimental MSDs of the depleted hydrated sample seem systematically higher than those of the α-La ca sample, which is hardly visible within the statistics in the simulation results at 280 K and below 1 ns. It indicates slightly enhanced dynamics for the depleted sample in such conditions, which could be expected as calcium has a stabilizing effect [START_REF] Permyakov | α-Lactalbumin: structure and function[END_REF]. The higher mobility becomes visible only in the simulations at higher temperatures and longer timescales, as the variations in the sample are certainly small. An interesting point is the difference between the models. For the simulations, the PK model mainly evaluates a very slightly smaller MSD values whereas for the experiments they were larger. One has to note that each spectrometer has not only its specific time resolution, but also a characteristic q-range. Both dimensions are important and are related. Therefore, the evaluated results do not only depend on the time resolution, but also on the accessible spatial domain, which permits us to see various behaviors of the samples.

Discussion and Conclusion

MD simulations are a very powerful tool to understand, in more detail, the dynamics of individual atoms that are measured for a sample in a neutron scattering experiment, as both techniques give access to comparable temporal and spatial scales. Unlike the common simulations run in solution, comparison to elastic incoherent neutron scattering (EINS) measurements, frequently done with hydrated powders, has required the development of approaches to simulate hydrated powders [START_REF] Yi | Derivation of Mean-Square Displacements for Protein Dynamics from Elastic Incoherent Neutron Scattering[END_REF][START_REF] Tarek | Environmental dependence of the dynamics of protein hydration water[END_REF] by adapting the setup accordingly.

A direct comparison of neutron data and simulated signals is not always trivial as the absolute values depend significantly on, one hand, data corrections and normalization, and on the other, on the accuracy of force fields and starting structures. It is also common to find that simulations cannot reproduce results extracted from neutron scattering data quantitatively (see [START_REF] Bachir Aoun | Direct comparison of elastic incoherent neutron scattering experiments with molecular dynamics simulations of DMPC phase transitions[END_REF] or Fig. 5.6), hence, the decision to compare MSDs by extracting them in a very similar way from both experiment and simulation. In addition, we checked that different approaches to calculate the static MSDs from the simulated trajectories gave identical results.

The order of magnitude of the values of experimental MSDs are well reproduced by the simulations, with MSDs of hydrated samples being larger than those of the dry samples. The results indicate that the models describing the simulated EISF (obtained from the DISF) underestimate the simulated directly calculated MSDs (see Fig. 5.4). One might therefore speculate that no model is able to take the whole dynamics into account and that effects due to the limited space and time windows are not negligible. For the hydrated protein, the differences are not as large as for the dry protein. In addition, the difference between the models is not negligible, but the trends are always the same and in agreement with the direct MSDs; i.e., all curves obtained through the different models are mainly parallel. Interestingly, the hierarchy between the models is always the same for an identical instrument resolution (with the exception of SPHERES in a dry environment). One can therefore conclude that the GA gives equally good results as the other models, since the absolute values of the MSDs are unknown.

In comparison to the experimental data, the simulation cannot provide reliable quantitative results (see Fig. 5.6). In general, the experimental MSDs of the PK model are higher than those from the GA model, whereas for the simulations this trend is inverted in most cases. Here, it has to be stressed again, this behavior is highly dependent on the chosen q-range and thus no definitive trend can be concluded. The experimental curves show larger differences and in particular the GA model gives MSDs which are more strikingly different from the MSDs obtained through the PK model. Nevertheless, none of these results favors any one model over another, as the statistics are probably not good enough to discriminate small effects, eventually due to the different q-ranges used.

As shown by Fig. 5.5 the distribution of the MSDs can be mainly described by two different peaks which are independent of hydration. The second peak is most visible above 1 ns, whereas below 60 ps it is not well distinguished. It is mainly the H atoms of the methyl groups (not shown here) that are contributing to this peak, which is in accordance to the findings of Yi et al. [START_REF] Yi | Derivation of Mean-Square Displacements for Protein Dynamics from Elastic Incoherent Neutron Scattering[END_REF]. Methyl group rotations indeed contribute to the elastic neutron spectra and the findings here support that they are a major contributor to heterogeneity originating from these motions, which becomes more visible at longer timescales. Yi et al. [START_REF] Yi | Derivation of Mean-Square Displacements for Protein Dynamics from Elastic Incoherent Neutron Scattering[END_REF] simulated the camphor-bound cytochrome P450 at h = 0.4 in a way comparable to the simulations here. They also showed that this peak is more dominant at higher temperatures. Furthermore, the second peak at larger amplitudes is also more pronounced at 1 ns. At 100 ps it is closer to the first peak and much broader. In addition, Tokuhisa et al. [START_REF] Tokuhisa | Non-Gaussian behavior of elastic incoherent neutron scattering profiles of proteins studied by molecular dynamics simulation[END_REF] simulated staphylococcal nuclease (SNase) in a water box at 300 K and also found two distinguishable peaks. The time was not documented but the evaluated simulation time was 1 ns, indicating that the investigated time window was likely smaller than 100 ps.

Overall this leads to the conclusion that the two models give reasonable results in comparison to the direct MSDs from the MD simulations. For a precise data set, the differences between the models are not significant concerning the trends, but the quantitative values are, depending on the evaluated q-range. The PK model gives further insight into the standard deviation of the MSD, but with respect to the MSD it does not give more accurate results. Furthermore, it is also important to state again that in contrast to the experimental data, the simulated EISF was not normalized to the lowest temperature data due to the lack of such simulation data, which could also partly explain the quantitative differences. Doing that, one would more consistently treat experimental and simulated data and eliminate more uncertainties, which might arise.

Chapter 6

Role of low-frequency vibrational dynamics of protein hydration water for ligand binding

Based on a paper in preparation:

Role of low-frequency vibrational dynamics of protein hydration water for ligand binding.

Daniele Di Bari, Judith Peters, Jacques Ollivier, Andrea Orecchini, Caterina Petrillo, Fabio Sterpone and Alessandro Paciaroni (to be submitted to Physical Review Letter)

We studied the low-frequency vibrational behavior of the hydration water of the maltose binding protein (MBP), a prototypical biomolecule where a large scale hingebending conformational change from the apo/open to the holo/closed state, the socalled Venus-flytrap mechanism, is key for ligand binding. By using inelastic neutron scattering spectroscopy we show that, upon complexation with maltose, not only the MBP, but also its hydration water displays significant low-frequency vibrational changes. The character of the MBP appears substantially softer in the apo than in the holo state, while the opposite is true for its hydration water. Normal mode analysis supports the experimental results and allows to identify their microscopic origin. The alternative energetic vibrational contributions of protein and protein hydration water support the view of the induced fit mechanism as underlying the MBP conformational switching.

Introduction

An accurate knowledge of the dynamics of proteins is necessary to deeply understand how they perform their biological activity [START_REF] Henzler-Wildman | Dynamic personalities of proteins[END_REF], a question that is fundamental to engineer proteins for specific functions and design next-generation therapeutics. Apart from the dominant contribution of large conformational changes and domain move-ments, the degree of protein flexibility which is relevant to a given biological process comes also from more subtle dynamical mechanisms [START_REF] Henzler-Wildman | Dynamic personalities of proteins[END_REF]. These motions include low-frequency (< 20 meV; < 5 THz) vibrational modes, the spectroscopic signatures of which have been measured in different functional states of proteins by inelastic neutron scattering (INS) [START_REF] Balog | Direct determination of vibrational density of states change on ligand binding to a protein[END_REF][START_REF] Niessen | Moving in the Right Direction: Protein Vibrational Steering Function[END_REF], optical Kerr effect spectroscopy [START_REF] Turton | Terahertz underdamped vibrational motion governs protein-ligand binding in solution[END_REF] and anisotropic terahertz microscopy [START_REF] Niessen | Moving in the Right Direction: Protein Vibrational Steering Function[END_REF][START_REF] Niessen | Protein and RNA dynamical fingerprinting[END_REF]. In the case of complexation, the protein vibrational density of states may either undergo softening [START_REF] Balog | Direct determination of vibrational density of states change on ligand binding to a protein[END_REF] or stiffening [START_REF] Niessen | Moving in the Right Direction: Protein Vibrational Steering Function[END_REF], with significant consequences in the increase or reduction of the complex flexibility. Quite remarkably, these changes affect, mainly through the vibrational entropy term, the stability of the protein-ligand complex and the binding affinity [START_REF] Balog | Direct determination of vibrational density of states change on ligand binding to a protein[END_REF][START_REF] Niessen | Moving in the Right Direction: Protein Vibrational Steering Function[END_REF]. On the other hand, the flexibility of a protein depends also on the dynamics of the water molecules interacting with its surface, the so-called hydration water (HW) [START_REF] Paul W Fenimore | Bulksolvent and hydration-shell fluctuations, similar to α-and β -fluctuations in glasses, control protein motions and functions[END_REF]. On the functional point of view, there exists already evidence that the HW diffusive dynamics plays a significant role in assisting enzyme-substrate interactions [START_REF] Grossman | Correlated structural kinetics and retarded solvent dynamics at the metalloprotease active site[END_REF][START_REF] Dielmann-Gessner | Enzymatic turnover of macromolecules generates long-lasting protein-water-coupled motions beyond reaction steady state[END_REF]. On the other hand, the role of the HW vibrational dynamics in protein-ligand recognition processes has been so far largely neglected. On these grounds we studied by INS spectroscopy the low-frequency vibrational behavior of the HW of the maltose binding protein (MBP), a prototypical member of a periplasmic-binding protein (PBP) superfamily [START_REF] Higgins | ABC TRANSPORTERS -FROM MICROORGANISMS TO MAN[END_REF]. Structurally, PBP proteins share a two-domain architecture (N-terminal domain, NTD, and C-terminal domain, CTD) with a central inter-domain cleft where the ligand is trapped following a large scale hinge-bending conformational change from the apo/open to the holo/closed state, the so-called Venus-flytrap mechanism. This conformational change is key for cellular metabolism [START_REF] Fa Quiocho | Atomic structure and specificity of bacterial periplasmic receptors for active transport and chemotaxis: Variation of common themes[END_REF], drug design [START_REF] Felder | The venus flytrap of periplasmic binding proteins: An ancient protein module present in multiple drug receptors[END_REF] and biosensor development [START_REF] Ribeiro | Converting a Periplasmic Binding Protein into a Synthetic Biosensing Switch through Domain Insertion[END_REF]. It is widely accepted that ligand recognition by MBP proceeds through an induced fit (IF) mechanism, by which the ligand binds the open state and prompts a transition to the closed state [START_REF] Lakshmi | Structure dictates the mechanism of ligand recognition in the histidine and maltose binding proteins[END_REF][START_REF] Marco Van Den Noort | Stability of ligandinduced protein conformation influences affinity in maltose-binding protein[END_REF], however the energetics behind this conformational change is still unclear.

Here we investigate the role played by protein HW vibrational changes upon complexation with maltose. We show that not only the very protein, but also its HW displays significant low-frequency vibrational changes, with the character of such changes being opposite in the MBP compared to its HW in terms of softening. The microscopic origin of this effect is dissected by exploiting normal mode analysis (NMA). Biomolecule and solvent have distinct and opposite contributions to the complexation vibrational free energy change.

Methods

Inelastic Neutron Scattering

INS experiments were done on the time-of-flight spectrometer IN5 (Institut Laue-Langevin, Grenoble) with an incident wavelength of 5 Å, (energy resolution ⇠ 0.09 meV) on apo-MBP (uncomplexed) and holo-MBP (complexed with maltose) hydrated powders. The most significant contribution to the measured spectra comes from the incoherent signal of the hydrogen atoms. In the case of the samples hydrated with heavy water, i.e. apo-and holo-MBP at h=0.40 g D 2 O/dry protein, the incoherent signal arises mainly from the protein non-exchangeable hydrogen atoms [START_REF] Varley | Neutron scattering lengths and cross sections[END_REF] that are copiously and almost uniformly distributed throughout the whole protein, thus allowing for a complete sampling of its molecular vibrational motions, within the time and spatial windows defined by the resolutions and ranges of the experimental energy transfer E and wave-vector transfer q. On the other hand, for apo-and holo-MBP hydrated with normal water, h=0.36 g H 2 O/dry protein, also the hydrogen atoms belonging to the HW give a substantial contribution to the measured spectra. The hydration degree h has been chosen to ensure that the average protein was in the presence of approximately the same number of H 2 O or D 2 O molecules in each sample corresponding to 1 -2 layers at the surface. The transmission of the samples, whose mass amounted to about 40 mg, ranged from 0.93 to 0.96. Multiple scattering and multiphonon contributions were estimated and considered negligible. The samples, placed in aluminum standard slab cells, were oriented at 135 . A temperature of 150 K was chosen, in such a way that anharmonic effects can be disregarded. Before any data analysis, the raw spectra were corrected for empty cell contribution and self-absorption, and normalized to a vanadium standard.

Normal Mode Analysis

All the simulations and the normal mode calculations were performed with GRO-MACS 2019.4 [START_REF] James | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF] with the CHARMM36m [START_REF] Yu | Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations[END_REF] force-field . The crystallographic coordinates of apo-MBP and holo-MBP were taken from the Protein Data Bank, entry 1OMP and 1ANF, respectively. The VMD software [START_REF] Humphrey | VMD: visual molecular dynamics[END_REF] and the python module parmed [START_REF] Michael R Shirts | Lessons learned from comparing molecular dynamics engines on the SAMPL5 dataset[END_REF] have been used to add the positions of the missing atoms, and generate the topology file for the ligand. Positions the crystal waters were not removed. Following the approach developed by Tarek and Tobias [START_REF] Tarek | The dynamics of protein hydration water: a quantitative comparison of molecular dynamics simulations and neutronscattering experiments[END_REF], using the PACKMOL software [START_REF] Martínez | PACKMOL: a package for building initial configurations for molecular dynamics simulations[END_REF], we created two boxes where we placed, at random, eight copies of the MBP proteins in the apo and holo forms, respectively, to reproduce the neutron experiments performed on an amorphous powder. Then, we hydrated the two systems with TIP3P H 2 O molecules to a hydration level of 36% (i.e. 36 g of H 2 O per 100 g of protein), and we added 64 Na ions in each box to neutralize the two systems.

To obtain the best staring configuration for the NMA, after a short energy minimization to eliminate unfavorable contacts, we performed a series MD simulations to relax the obtained structures. The particle mesh Ewald (PME) method was employed to compute the long-range electrostatic interactions, and a cutoff distance of 1.2 nm was used for calculation of the short-range electrostatic and van der Waals interactions. For the minimization, the steepest descent algorithm was used to a tolerance of 1000 kJ/mol/nm. The MD simulations were carried out with an integration time step of 1 fs. The systems were first heated to 300 K with a temperature increment of 10 K after each 50 steps. At each step, a temperature equilibration (NVT) was followed by a constant pressure and temperature MD run (NPT) to allow the proteins to interact with their neighbors and periodic images, leading to a contraction of the boxes. After a longer NPT equilibration at 300 K of 100ns, the systems were cooled down to 150 K with a reduction of 10 K every 50 ps. A further long NPT equilibration was performed until the system reached a stationary volume. To keep constant the required temperature in the NVT and the NPT runs, we used the V-Rescale thermostat, coupled with the Parrinello-Rahman barostat to maintain the pressure at 1 bar.

At the end of the last NPT run, the coordinates of each of the eight proteins, both for the apo and holo systems, were taken separately with their closest 791 H 2 O molecules (36% hydration level). Each hydrated protein was energy minimized with the steepest descent algorithm [cit] followed by the low-memory Broyden-Fletcher-Goldfarb-Shanno method until the maximum force in the system was smaller than 10 4 kJ/mol/nm.

Discussion

A direct view of the way the low-frequency modes of MBP and its hydration water are affected by ligand binding is provided by Fig. 6.2, where the dynamic structure factor S(q, E) as measured by inelastic neutron scattering is reported. Panel (b) shows that the dynamic structure factor of holo-MBP displays an excess of inelastic signal compared to apo-MBP, This behaviour is in line with that reported for dihydrofolate reductase protein [START_REF] Balog | Direct determination of vibrational density of states change on ligand binding to a protein[END_REF], but deviates from the low-frequency trend found for lysozyme [START_REF] Turton | Terahertz underdamped vibrational motion governs protein-ligand binding in solution[END_REF][START_REF] Niessen | Moving in the Right Direction: Protein Vibrational Steering Function[END_REF] and aspartate aminotransferase [START_REF] Dajnowicz | Substrate Binding Stiffens Aspartate Aminotransferase by Altering the Enzyme Picosecond Vibrational Dynamics[END_REF] where ligand binding is related to a rigidified vibrational dynamics. Due to the close interaction of water molecules with the protein surface, one may ask whether there will be a vibrational change upon complexation also in the protein HW dynamic structure factor. The hydration water spectra have been obtained by properly subtracting the contribution of apo-MBP+D 2 O (holo-MBP+D 2 O) from the signal of MBP-APO-H 2 O (MBP-HOLO-H 2 O), to obtain the HW S(q, E) for the free (complexed) protein (See the SI for details). Panel (c) of Fig. 6.2 shows that the apo-MBP HW displays an excess of inelastic signal compared to holo-MBP HW, a behavior opposite to that of the protein. Indeed, to quantitatively describe the vibrational changes occurring upon complexation to both the protein and its HW, it is convenient to calculate the corresponding vibrational density of states g(E), which is directly related to the dynamic structure factor by g(E)=lim q!0 6E }q 2 (exp E K B T 1)S(q, E). Actually, the g(E) is mainly the proton weighted vibrational density of states, as it derives from the average over the strong incoherent signal from all the hydrogen atoms. The vibrational modes of the g(E) shown in panel (c) of Fig. 6.2 arise mostly from the collective displacements of both protein groups of atoms and water molecules, while in panel (b) and panel (c) the protein and HW contributions are singled out respectively. The g(E)s have been normalized to their absolute values using the procedure described in Ref. [START_REF] Balog | Direct determination of vibrational density of states change on ligand binding to a protein[END_REF] (See the SI for details). A vibrational difference beyond the errorbars between the free and the complexed protein is visible in the low-energy range from 2 meV to 4 meV, with a higher number of low-frequency large-amplitude modes appearing in the holo state. This excess of modes provides complexed MBP with additional flexibility, a feature also known as softening. On the other hand, the protein HW displays a reverse behaviour, with the free system showing an excess of low-frequency modes compared to the case in the presence of maltose. To explain the origin of the antithetical trend shown by the protein and HW g(E) we examined the structural changes occurring to MBP during the apo/open to holo/closed transition. In particular, a key role for this transition is played by protein hinge region (residues K170-D180), which includes a short α-helix and a two-stranded β -sheets NTD [START_REF] Millet | The energetic cost of domain reorientation in maltose-binding protein as studied by NMR and fluorescence spectroscopy[END_REF]. The greater flexibility of this region after complexation, also revealed by the increased average B-factors [START_REF] Patrick | Insights into the conformational equilibria of maltose-binding protein by analysis of high affinity mutants[END_REF], is related to the excess of low-energy excess of modes displayed by holo-MBP. Indeed, numerical g(E) non only reproduces the trend observed by experiments (Fig. 6.3), but also confirms that the low-frequency vibrational density of states of the atoms in the hinge region is larger in the holo than in the apo conformation and is predominant over the remaining protein contribution (see Fig. 6.4). On the other hand, Fig. 6.4 shows that the apo-MBP HW is more ordered than the holo-MBP, and this means that it has a lower density [START_REF] Kt Wikfeldt | Spatially inhomogeneous bimodal inherent structure of simulated liquid water[END_REF].

On the thermodynamic point of view, the features in the g(E) reflect directly into the changes of the vibrational free energy A vib of both the protein and its hydration water on passing from the apo/open to the holo/closed conformation. In the harmonic approximation one can easily estimate such changes through the relationship [START_REF] Maradudin | Theory of lattice dynamics in the harmonic approximation[END_REF]:

∆A α vib = A α vib,holo A α vib,apo = = K B T Z ln  2 sinh ✓ E 2K B T ◆ g α holo (E) g α apo (E) dE (6.1)
where α = MBP or HW . As for the protein, it turns out that ∆A MBP vib =-(20 ± 3) kJ/mol, i.e. the vibrational change provides a favorable contribution for complexation. On the other hand, and quite surprisingly, the contribution relative to the MBP HW ∆A HW vib = (14 ± 5) kJ/mol promotes the stabilization of the open state. These results are in agreement with the fact that just the MBP intrinsic dynamics is insufficient for a 'selection' mechanism [START_REF] Kim | A single-molecule dissection of ligand binding to a protein with intrinsic dynamics[END_REF], thus supporting the IF model for this enzyme. The emerging picture is one where substrate-mediated interactions are necessary to induce the protein open-to-closed conformational transition and overcome the energy barrier separating the MBP apo/open and holo/closed conformational state [START_REF] Lakshmi | Structure dictates the mechanism of ligand recognition in the histidine and maltose binding proteins[END_REF] to which the protein HW vibrational dynamics contributes significantly through the ∆A HW vib term. The holo/closed state is further stabilized by the excess of low-frequency protein modes via the contribution ∆A MBP vib value, mainly of entropic nature. Quite remarkably, both ∆A HW vib and ∆A MBP vib are of the same order of the total complexation free energy, which amounts to 35 ± 1 kJ/mol [START_REF] Millet | The energetic cost of domain reorientation in maltose-binding protein as studied by NMR and fluorescence spectroscopy[END_REF], thus indicating that protein and HW vibrational degrees of freedom represent a crucial part in the energetics of the ligand binding processes. where u is the displacement of the atom at the l-site of the lattice -i.e. u is the fluctuation of the position r(t) of an atom with respect to its average position hr(t)i = l. Then, we can consider the relation between the Debye-Waller factor and the atom displacement (see eq. (4.38) -Lovesey [START_REF] Stephen W Lovesey | Theory of neutron scattering from condensed matter[END_REF]):

2W (κ)= ⌦ (κ • u) 2 ↵ = κ 2 ⌦ u 2 κ ↵ (6.4)
where ⌦ u 2 κ ↵ is the mean value of the projection of the atomic fluctuation into the direction of κ -i.e. if we define the direction of κ as n κ , u κ is defined as:

u κ = n κ • u
Therefore, combining eq. 6.2 and eq. 6.4 we get:

⌦ u 2 κ ↵ = 3 h 2 M κ 2 Z ω m 0 dω g(ω) ω • coth ✓ h ω 2 k b T ◆ • n |κ • σ | 2 o avg (6.5)
Now, in the case of cubic symmetry (or disordered materials), we have: Now, if we take the limit of t 0 ! ∞, we have that r(t + t 0 ) and r(t 0 ) are uncorrelated: lim t 0 !∞ hr(t + t 0 ) • r(t 0 )i t 0 = hr(t + t 0 )i t 0 • hr(t 0 )i t 0 = hri 2 and consequently, from eq. 6.13, we obtain:

n |κ • σ | 2 o avg = 1
MSD(t ! ∞)=2 ⇣ ⌦ r 2 ↵ hri 2 ⌘ ⌘ 2 MSPF (6.14)
Therefore, combining eq. 6.9, 6.10, and 6.14 we obtain:

2W (κ)= 1 6 κ 2 MSD(t ! ∞) (6.15) MSD(t ! ∞)= 3 h M Z ω m 0 dω g(ω) ω • coth ✓ h ω 2 k b T ◆ (6.16)
In conclusion, if we want to express eq. 6.10 and 6.16 in function of E instead of ω, we can make the substitutions:

hω ! E dω g(ω) ! dE g(E)

such that 1 :

MSPF = 3 h2 2 M Z E m 0 dE g(E) E • coth ✓ E 2 k b T ◆ (6.17) MSD(t ! ∞)= 3 h2 M Z E m 0 dE g(E) E • coth ✓ E 2 k b T ◆ (6.18)

Normalization on the vDOS

The vDOS obtained from the incoherent dynamic structure factor (IDSF) is on relative scales. To re-normalize the vDOS we can use the ratio between the experimentally measured mean square MSD, which is on absolute scale, and the one calculated from the vDOS on relative scale. Practically, this means that the vDOS in absolute scale, g abs (E), can be calculated as:

g abs (E)= MSD exp MSD calc • g(E) (6.19)
where MSD exp is the value of the MSD measured from the elastic peack at low q (i.e. from S(q, E = 0)). In the gaussian approximation:

S(q, E = 0)=S 0 • e 1 6 q 2 • MSD exp ! MSD exp = 1 6 • ∂ ln S(q, E = 0) ∂ q 2 (6.20)

However, the our IDSF data were collected at IN5 that is not optimized for this type of elastic measurements. Therefore, we decided to take, for the value of the MSD exp of the apo-MBP, the one measured by Wood et al. [START_REF] Wood | Coincidence of dynamical transitions in a soluble protein and its hydration water: direct measurements by neutron scattering and md simulations[END_REF], and in order to re-normalize properly also the vDOS of the holo-MBP, we considered that the differences between in the vDOS of apo-MBP Apo and holo-MBP actually are due to a shift of the vibrational modes due to the complexation [START_REF] Balog | Direct determination of vibrational density of states change on ligand binding to a protein[END_REF]. This means that, if there is a softening of the complexed MBP at low energy (vDOS holo > vDOS apo ), there must be a stiffening at higher frequencies (vDOS holo < vDOS apo ) that compensates for the softening. Consequently, if we assume that this compensation takes place before of a certain energy E, then the integrals of the vDOS between 0 and E of the Apo and the Holo systems should be equals and we can use their ratio to re-scale the vDOS. where we assumed that the compensation take place before 15 meV, i.e. that E = 15 meV. Finally, with the vDOS properly normalized, we were able to obtain the vDOS of the HW, for both the apo and the holo systems, simply subtracting the vDOS of measured for the MBP from one measured from MBP+H 2 O.

vDOS from NMA

The vibrational density of states of a system is the frequency distribution of the vibrational modes, obtained by solving the eigenvalue problem for the mass-weighted Hessian matrix [START_REF] Balog | Vibrational softening of a protein on ligand binding[END_REF]

g α (E)= 1 ∆E 3N ∑ i=1 |e i (α)| 2 δ (E E i ) (6.23)
where e i (α) is a three-dimensional vector with the x, y, and z components of the massweighted displacement of atom α within the normal mode i (i.e. it is the eigen-vector obtained with the NMA), ∆E is the width of the sampling interval and δ (E E i ) such that, it is equal to 1 when ∆E 2  E E i < ∆E 2 , otherwise it is 0.

Domain Opening Angle (DOA)

The MBP DOA is defined as the angle between the Centers of Mass (CoM) of the C atoms of the Nt domain (i.e. the N-terminal residues 1-108, 263-311), the CoM of the C of central t domain (i.e. the -sheet hinge residues: 109-112, 259-262) and the COM of the C of the Ct domain (i.e. the C-terminal residues 113-258, 320-370).

Tetrahedral Order Parameter

The local tetrahedrality is a many-body property and can be defined (Chau & Hardwick, 1998) by:

q i = 1 3 8 • 3 ∑ j=1 4 ∑ k= j+1 ✓ cos θ jik + 1 3 ◆ 2 (6.26)
where the sum is over all angles θ jik formed around a reference molecule i by its four nearest potentially hydrogen-bonding neighbors, although no account is taken of whether the neighbors are actually hydrogen bonded to molecule i. A larger value of q i indicates a greater local tetrahedrality. In the case of water molecules close to the protein, nearest neighbors are taken to include potential hydrogen-bonding moieties of the protein, such as threonine hydroxyl groups or the amide bonds of the protein backbone. Close to the protein surface, there will undoubtedly be water molecules in a constricted environment which cannot form tetrahedral structures.

Local structure index (LSI)

Let r i, j be the radial distance between the oxygen of the molecule i and the oxygen of molecule j ordered such that r 1 < r 2 <...<r i, j <...<r i,n(i) < r i,n(i)+1 , where n(i) is chosen so that r i,n(i) < 3.7 Å < r i,n(i)+1 . Then, the LSI parameter is defined by as (Shiratani & Sasai, 1996 and 1998):

LSI i = 1 n(i) • n(i) ∑ j=1 ⇥ ∆ i, j ∆ i ⇤ (6.27)
where ∆ i, j = r i, j+1 r i, j , ∆ i is the average value of ∆ i, j over all the n(i)+1 molecules.

The LSI aims at measuring the extent of the gap between the first and the second hydration shells surrounding a water molecule measuring the inhomogeneity in the radial distribution within the sphere of radius 3.7 Å. A high value of LSI i implies that the molecule i at time t is characterized by a tetrahedral local order and a low-local density, while on the contrary, values of LSI i ⇡ 0 indicate a molecule with defective tetrahedral order and high-local density [START_REF] Gustavo | Evidence of a two-state picture for supercooled water and its connections with glassy dynamics[END_REF].

Conclusions and discussion

The work presented in this thesis aims at shedding some microscopic insights into thermal stability of bacterial cells (denaturation and cell growth). In particular, we found that the thermal death of E. coli is signaled by a distinctive behavior of the proteome short-time dynamics: a strong decrease of the protein global diffusion coefficient starts just below T CD , from D G (320 K)= 1.5 Å2 /ns down to D G (350 K)= 0.5 Å2 /ns. The results from MD simulations show that this dynamic slow-down is due to the unfolding of a part of the proteome. This finding is consistent with previous experimental work where an analogous trend for the diffusive dynamics has been observed in concentrated protein solutions across the melting temperature [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF][START_REF] Grimaldo | Hierarchical molecular dynamics of bovine serum albumin in concentrated aqueous solution below and above thermal denaturation[END_REF][START_REF] Matsarskaia | Evolution of the structure and dynamics of bovine serum albumin induced by thermal denaturation[END_REF]. Therefore, protein unfolding dominates the observed temperature dependence, and in both cases the dynamic slow-down is irreversible.

In this study, we take the analysis a crucial step forward by directly relating the diffusion coefficient to the amount of unfolded proteins in the system. The theoretical description of protein solutions containing a varying ratio of folded and unfolded proteins is challenging. These systems neither fit into the standard picture of solutions formed by globular proteins, usually modeled as rigid colloidal particles [START_REF] Stradner | Potential and limits of a colloid approach to protein solutions[END_REF], nor can they be accurately described by means of simplified polymer models [START_REF] Ralph | Structure and linear viscoelasticity of flexible polymer solutions: comparison of polyelectrolyte and neutral polymer solutions[END_REF]. In particular, limitations of the colloidal model were demonstrated for highly concentrated protein solutions involving changes in protein conformation [START_REF] Prasad S Sarangapani | Critical examination of the colloidal particle model of globular proteins[END_REF], and the importance of protein-protein interactions was stressed [START_REF] Douglas Godfrin | Generalized phase behavior of cluster formation in colloidal dispersions with competing interactions[END_REF][START_REF] Zhang | Recent progresses of understanding the viscosity of concentrated protein solutions[END_REF].

Our simulations showed that the presence of minor amounts of unfolded proteins causes a substantial slow-down in the protein global diffusion. As we demonstrated, this drop was not only due to the slower diffusion of the fraction of unfolded proteins, but also to a twofold slow-down of the remaining folded proteins as a consequence of the enhanced interactions with their unfolded counterparts. Thus, unfolded proteins form a sticky macromolecular network to which folded proteins associate. This resembles the behavior recently observed in biomolecular condensates, where the interactions of folded lysozyme proteins with a macromolecular network formed by pentameric constructs of SH3 domains, and containing disordered linkers, strongly affected the condensate viscoelastic properties [START_REF] Ghosh | Shear relaxation governs fusion dynamics of biomolecular condensates[END_REF].

From the combination of QENS experiments and MD simulations, we estimated the amount of unfolded proteins in the cytoplasm at different temperatures. In the last years, there have been several attempts to connect the thermal death of bacteria to a critical amount of unfolding proteins. The aim was to understand if the death results from a collective unfolding of cell proteins [START_REF] Ghosh | Cellular proteomes have broad distributions of protein stability[END_REF][START_REF] Dill | Physical limits of cells and proteomes[END_REF] or if it is caused by the denaturation of a subset of proteins controlling key biological functions [START_REF] Roger L Chang | Structural systems biology evaluation of metabolic thermotolerance in Escherichia coli[END_REF][START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF][START_REF] Mateus | Thermal Proteome Profiling in Bacteria: Probing Protein State in vivo[END_REF][START_REF] Jarzab | Meltome atlas-thermal proteome stability across the tree of life[END_REF]. Here, we found that a few degrees above the cell-death temperature only a small fraction of proteins, less than 15%, are unfolded (see also Table 4.10). This result supports the hypothesis firstly put forward by Leuenberger et al. [START_REF] Leuenberger | Cell-wide analysis of protein thermal unfolding reveals determinants of thermostability[END_REF] that there is no catastrophic denaturation of the proteome, but instead only an unfolding of a subset of proteins. It is important to stress that there is no unique definition of T CD which, depending on the growth conditions of bacteria and their environment, can vary by several Celsius degrees. Owing to the predicted slow increase with temperature of the unfolded protein fraction (see Fig. 4.28), this uncertainty does not affect our conclusions concerning the minor amount of unfolded proteins that are present in the cytoplasm at the cell death. Apart from the uncertainty in T CD , the quantitative estimate of the unfolded fraction may also be affected by some limitations of our computational model in terms of molecular composition, which focuses exclusively on proteins as most prevalent type of macromolecules in the cytoplasm and which, moreover, is biased toward structurally well-resolved folded proteins. In addition, the model does not consider thermal adaptations of the proteome, such as evolving populations of heat-shock proteins and chaperons. However, even though the cells in the experimental samples still have a basal metabolism, they will have a reduced capability to tune the proteome composition due to lack of nutrients.

The destabilization of the E. coli bacteria starts already at temperatures below the T CD . Temperatures near 315K represent already a stress condition for the bacteriathey resist to the increase of the environmental temperature with several active mechanisms, such as the change in the global protein population by increasing the number of molecular chaperones to maintain a properly folded proteome [START_REF] Chen | Thermosensitivity of growth is determined by chaperonemediated proteome reallocation[END_REF] and the variation of internal viscosity by regulating the synthesis of glycogen and trehalose [START_REF] Laura B Persson | Cellular control of viscosity counters changes in temperature and energy availability[END_REF]. With a contribution of such mechanisms of thermal adaptation, the average dynamical state of the cytoplasm observed at T CD is still similar to the dynamical state of the cytoplasm near the temperature of optimal growth rate. This leads us to argue that the average dynamics should not play a predominant role in thermal death. However, as our simulations have shown, a low amount of unfolded proteins can trigger an important slow-down in the diffusion of the surrounding macromolecules. Therefore, the increase of local viscosity and the associated dramatically reduced protein diffusion caused by unfolding, may threaten the viability of the cell by affecting localized physiological processes. A pertinent example at molecular scale is the dynamics and substrate channeling in enzymatic assembling [?], but also at larger scale the viscoelastic response of the cytoplasm associated to organelles localization [START_REF] Xie | Contribution of cytoplasm viscoelastic properties to mitotic spindle positioning[END_REF]. Morever, we have shown that a very good reproduction of E. coli growth-rate can be obtained when combining together the fraction of unfolded proteins with the temperature-dependent diffusion within a simple reaction-diffusion model. This supports the idea that at least a part of the cell metabolism is modulated by diffusion. Finally, our findings are consistent with the idea that intrinsically disordered proteins may induce as well a mobility slow-down on the local environment, as it is probed in membrane-less organelles [START_REF] Taylor | Quantifying Dynamics in Phase-Separated Condensates Using Fluorescence Recovery after Photobleaching[END_REF].

The approach hereby presented can be extended -and possibly complemented by single-molecule techniques-to investigate the relationship between the dynamics and the proteome unfolding in extremophiles resisting either to cold or hot environments. Further, attention could be enlarged to the peculiarity of the proteins' dynamical response to stress in the functioning of specific networks of interaction like in the unfolding protein response cascade.

In this thesis we presented also two other studies: "Differences between Ca 2+ reach and depleted α-La investigated by MD simulations and NS experiments", and "Role of low-frequency vibrational dynamics of protein hydration water for ligand binding"dealing with the effects of protein-ligand complexation on protein dynamics, and in which neutron scattering techniques and molecular dynamics simulations are always coupled.

In the first one, we show how data from EINS and MD can be combined to probe the dynamics in both protein-depleted and protein calcium rich systems. Two models were exploited to extract the hydrogen atoms MSD from the EINS data, measured at three different levels of hydration and with resolutions to explore different time scales (picosecond to nanosecond, i.e. a time scale range accessible to simulations). In parallel, systems mimicking protein powders at different hydration levels were simulated at the atomistic level, and MSD estimated with different methods, allowing a robust comparison with neutron scattering data.

In the last study presented in this thesis, we describe the behavior of hydration water in a system containing MBP, a protein undergoing a large-scale conformational change (open-to-close) in ligand binding. In this case, INS experiments were performed to probe the low frequency vibration modes of protein and water in both states. Results point out that the protein and its hydration water have an opposite behavior in closing. An interesting thermodynamics scenario, which however deserves further studies, is proposed where the protein-ligand interactions allow the protein to overcome the opento-close free energy barrier, which is contributed by the hydration water (in a way still to be clarified).

Figure 1 . 1 :

 11 Figure 1.1: Common features of living cells. Every cell has a nucleus or nucleoid storing their DNA, a plasma membrane, and cytoplasm. Source: D. L. Nelson and M. M. Cox, "Lehninger Principles of Biochemistry" (7th edition, 2017) [3].

Figure 1 . 2 :

 12 Figure 1.2: Transmission electron microscopy (TEM) images of E. coli cells showing the damage of their membranes. A) incubated at 35 C, and heated to B) 50 C, and C) 65 C. The first row images are taken with x13,500 magnification; second row images are taken with x92,000 magnification. Source: B. Tonyali et al., "Evaluation of heating effects on the morphology and membrane structure of Escherichia coli using electron paramagnetic resonance spectroscopy" (Biophysical Chemistry 2019) [30].

Figure 1 . 3 :

 13 Figure 1.3: Hydrodynamic diameter of E. coli cells with heat treatment (35-70 C), measured by dynamic light scattering (DLS). The solid line represents the mean of individual particle size measurements and the dash lines represent standard deviations. Source: B. Tonyali et al., "Evaluation of heating effects on the morphology and membrane structure of Escherichia coli using electron paramagnetic resonance spectroscopy" (Biophysical Chemistry 2019) [30].

  .4.

Figure 1 . 4 :

 14 Figure 1.4: Proteome reallocation with change in temperature. The percentage is calculated with FoldME. Source: K. Chen et al., "Thermosensitivity of growth is determined by chaperonemediated proteome reallocation" (PNAS 2017) [43].

Figure 2 . 1 :

 21 Figure 2.1: Neutron-induced fission of the 235U nucleus. Source: [54].

Figure 2 . 2 :

 22 Figure 2.2: The flux per unit energy from a moderator at three temperatures.

Figure 2 . 3 :

 23 Figure 2.3: Schematic of the high-flux reactor at the Institut Laue-Langevin in Grenoble, France. Labels of the form 'Hn' indicate beam tubes and guides. Source: [54].

Figure 2 . 4 :

 24 Figure 2.4: Geometry for a scattering process. Source: [51].

Figure 2 . 5 :

 25 Figure 2.5: Illustration of the EINS, QENS and INS component of S inc (q, E). Source: [55].

Figure 3 . 1 :

 31 Figure 3.1: A schematic representation of phase space. The hexagonal cells represent state points (r N , p N ). In an ergodic system, all the trajectories in this figure represent different sections of a single long trajectory. Indeed, if the system is ergodic, the single long trajectory would eventually pass through (or arbitrarily near) all states. A substantial region of cyclical trajectories, and a barrier region leading to bottleneck, are shaded. Source: Allen and Tildesley, Computer Simulation of Liquids (1st edition, 1987) [77].

Figure 3 . 2 :

 32 Figure 3.2: Schematic representation of the bonded interaction terms contributing to the force field. Source: P. Gkeka and Z. Cournia, Molecular Dynamics simulations of lysozyme in water (2015) [79].

Figure 3 . 3 :

 33 Figure 3.3: Schematic representation of Lennard-Jones potential. The collision parameter, σ , is shown along with the well depth, ε. Source: Eni. Generalic, Lennard-Jones potential (Croatian-English Chemistry Dictionary & Glossar, 2017) [80].

Figure 3 . 4 :

 34 Figure 3.4: A two-dimensional periodic system. Molecules can enter and leave each box across each of the four edges. In a three-dimensional example, molecules would be free to cross any of the six cube faces. Source: Allen and Tildesley, Computer Simulation of Liquids (1st edition, 1987) [77].

Figure 3 . 5 :

 35 Figure 3.5: In PME, a charge (denoted by an empty circle with label q in the figure) is distributed over grid (here a mesh in two dimensions) points with weighting functions chosen according to the distance of the respective grid points to the location of the charge. Positioning all charges on a grid enables the use of the Fast Fourier Transform (FFT) to solve Poisson's equation for the electrostatic potential due to the charge distribution. In real applications, the grid is three-dimensional. Source: J. C. Phillips et al., Scalable Molecular Dynamics with NAMD (Journal of Computational Chemistry, 2005) [81].

Figure 3 . 6 :

 36 Figure 3.6: Simple example that shows the merit of a symplectic integrator: integration of a one-dimensional harmonic oscillator with an unit circle as exact trajectory. The trajectory of nonsymplectic method initially draw a circle, but after few steps starts to collapse toward the center while the symplectic method maintains a stable orbit, showing a superior long-time stability, even though its trajectory is deformed into an ellipse by a larger local error. Source: J. C. Phillips et al., Scalable Molecular Dynamics with NAMD (Journal of Computational Chemistry, 2005) [81].

Figure 3 . 7 :

 37 Figure 3.7: Schematic representation of the velocity Verlet algorithm. At each step, the stored variables are in grey boxes. Source: Allen and Tildesley, Computer Simulation of Liquids (1st edition, 1987) [77].

Figure 3 . 8 :

 38 Figure 3.8: Simplified example of the energy landscape (as a function of only two variables). Source: Voet D., Voet J.G. and Pratt C.W., "Fundamentals of Biochemistry" (5th edition 2016)[83].

Figure 3 . 9 :

 39 Figure 3.9: From the Mountains to the Valleys: a molecular dynamics fairy tale. Source: Theoretical and Computational Biophysics Group, "Computational Biophysics Workshop" (Boston, Dec. 5-9, 2004).

Figure 3 . 10 :

 310 Figure 3.10: Schematic representation of various computational approaches, used to cover different scales of length and time pertinent to different biophysical processes. These methods range from highly accurate, but computational demanding to the highly efficient but very low-detail continuum models. Figure reprinted from [84].

Figure 3 . 11 :

 311 Figure 3.11: All-atom versus coarse-grained energy landscape. The figure illustrates the effect of the smoothening of the energy landscape in a coarsegrained model as compared to an all-atom model. Figure reprinted from [84].

Figure 3 . 12 :

 312 Figure 3.12: The figure show an example of the OPEP model for the peptide Ala-Lys-Phe-Pro-Val in its zwitterion form to show the details of the backbone and the side-chains. Figure reprinted from [85].

Figure 3 . 13 :

 313 Figure3.13: Pictorial view of the LBMD multiscale scheme. Proteins are described at microscopic level, interact according to the OPEP coarsegrained force field and move in the continuum. The aqueous solvent is handled by the lattice Boltzmann method, so that fluid populations that reside on a Cartesian mesh are evolved in time and move to neighboring mesh points as connected by a set of discrete speeds (indicated as D3Q19)[START_REF] Chiricotto | Hydrodynamic effects on β -amyloid (16-22) peptide aggregation[END_REF].

Chapter 4 Characterization

 4 of the Dynamical State of the E. Coli Cytoplasm and the Effect of Cell Death Based on a paper in preparation: Short-time Diffusive Dynamics of Bacterial Proteome as a Proxy of the Cell Death. Daniele Di Bari, Stepan Timr, Marianne Guiral, Marie-Thérèse Giudici-Orticoni, Tilo Seydel, Christian Beck, Caterina Petrillo, Philippe Derreumaux, Simone Melchionna, Fabio Sterpone, Judith Peters and Alessandro Paciaroni (to be submitted to Science)

Figure 4 . 1 :

 41 Figure 4.1:Example of the QENS spectrum measured for the Vanadium at 275K and q = 0.94 Å 1 . The red line is the best fit of R(q, E) described by the eq. (4.3), meanwhile the dashed lines are the different components of R(q, E). The inset shows the Full Width Half Maximum (FWHM) of R(q, E) which represents a measure of the energy resolution of the instrument.

Figure 4 . 2 :

 42 Figure 4.2: PBS-D2O buffer. Resulting parameters of the fit with the eq. (4.4).

Figure 4 . 3 :

 43 Figure 4.3: Global motions (1 st step). Lorentzian widths γ G vs. q 2 for different temperatures. The dashed lines are fits with the jump-diffusion model described by the eq. (4.6).

Figure 4 . 4 :

 44 Figure 4.4: Local motions (1 st step). Lorentzian widths γ L vs. q 2 for different temperatures. The dashed lines are fits with the jump-diffusion model described by the eq. (4.6).

Figure 4 . 5 :

 45 Figure 4.5: Diffusion coefficients and residence times for the global and local motions obtained from the simultaneous fit of the E. coli data, assuming a priori the jump-diffusion model eq. (4.6). The solid blue line is a fit of D G with the gelation model described by the eq. (4.7)

  21 ± 0.06) C a 1 (0.020 ± 0.001) Å2 /(ns •C) b 1 (0.605 ± 0.005) Å2 /ns a 2 10 10 ± 0.001 Å2 /(ns •C) b 2 (0.48 ± 0.01) Å2 /ns φ (40.1 ± 0.1) %An example of the resulting parameters for the EISF, A 0 (q), at five selected temperatures is shown in Fig.4.6 (left).

Figure 4 . 6 :

 46 Figure 4.6: Left: EISF derived from the global fit of the E. coli data. Dashed lines are the best fit of the EISF with the model described by eq. (4.9). Right: Parameters resulting from the fit of the EISF. p L is the fraction of H-atoms appearing fixed on the accessible time scale of the instrument, s L is the fraction of H-atoms diffusing in a spherical volume with radius r L .

Figure 4 . 7 :

 47 Figure 4.7: EINS MSD fitted with eq. (4.7) describing the gelation of the system. The resulting parameters, T 0 and ∆T , are shown in the box on the upper left. These values are compared with the corresponding parameters obtained from QENS (box on the upper right.

  .

  8)were converted to the OPEP description (see section 3.3.2) by the OPEP File Generator available on http://opep.galaxy.ibpc.fr. Initial positions and orientations of the proteins in the simulation box were generated by the Packmol software[START_REF] Martínez | PACKMOL: a package for building initial configurations for molecular dynamics simulations[END_REF]. To make the spatial distributions of different protein species more uniform, each protein was treated separately in Packmol.

Figure 4 . 8 :

 48 Figure 4.8: Structures and counts of proteins included in the LBMD simulation.

Figure 4 . 9 :

 49 Figure 4.9: Dependence of the LBMD coupling coefficients γ on protein molecular weight. The γ values for six proteins (CspC, SOD1, Adk, TufA, MetE, and GapA) were optimized in order for the diffusion coefficient of the given isolated protein in an LBMD simulation to match the HYDROPRO [137] prediction. To determine the γ values for the remaining 29 protein species, the resulting dependence of 1/γ on molecular weight (MW) was fitted with a power-law dependence 1/γ = a • MW b + c, yielding the parameter values a = 35.76222504, b = 0.81676473, and c = 113.05584831 (for MW expressed in kDa and γ in fs 1 ).

Time step: 1 fsFFigure 4 . 10 :

 1410 Figure 4.10: Schematic representation of the procedure used to equilibrate atomistic subvolumes at different temperatures and in different folding states. (A) General overview of the procedure, (B-F) detailed description of its individual parts.

Figure 4 . 11 :

 411 Figure 4.11: Monitoring the unfolding protocol for each protein in the sub-volume with the concentration of 287.7 g/L. (A) Root-mean-square deviation, RMSD, of the proteins. (B) Modulus of radius of gyration. (C) Number of residues per secondary structure, as determined by the DSSP software: the graph in the upper right corner is a magnification of other main plot on the left to focus on the effects of the heat-shock at 1500 K.

Figure 4 . 12 :

 412 Figure 4.12: Examples of MSD curves in a folded-(blue) and an unfolded (red) system. The MSD curves were computed as a weighted average over all protein center-of-mass MSDs in 20 ns blocks of a production trajectory performed at T = 300 K for the 288 g/L sub-volume using the Amber99SB-disp force field.

Figure 4 .

 4 Figure 4.13:The η ij (t) running integrals (Eq. (4.14)) obtained for the 288 g/L CHARMM36m sub-volume at T = 300 K. The black curve shows the average of η ij (t) over all the NV T simulations, the gray stripe represents the standard error of the mean, and the blue curve is a fit using the analytical function(4.15).

Figure 4 .

 4 14, and the corresponding values of the parameters B w , D w , and ∆E w are listed in Table 4.8.

Figure 4 . 14 :

 414 Figure 4.14: Solvent viscosities determined from simulations. The solvent was modeled as a 150 mM aqueous KCl solution, and the respective water model (Amber99SB-disp water, or TIP3P) was used. Intermediate viscosity values (blue line) were interpolated using the expression (4.16).

Figure 4 . 15 :

 415 Figure 4.15: Viscosities calculated for the 288 g/L sub-volume using the Amber99SB-disp force field (left column) and the CHARMM36m force field (right column). (A) Folded system. (B) Unfolded system. The solid lines show an interpolation using Eq. (4.17) with bands expressing the uncertainty of this interpolation. (C) Extrapolation to different concentrations using Eq. (4.17).

Figure 4 . 16 :

 416 Figure 4.16: Comparison of D for the folded and the unfolded average protein at different temperatures for systems with different concentrations (force field: CHARMM36m).

Figure 4 .

 4 Figure 4.17: D of the average protein at different temperatures for systems with different concentrations (force field: CHARMM36m). The slowdown is reproduced using the smeared step function a QENS u to weight the prediction obtained from the simulation of the fully folded and fully unfolded systems.

Figure 4 .

 4 Figure 4.18: A: Probability distribution for the selection of a sub-volume of (17 nm) 3 with a concentration c from the larger system of (40 nm) 3 used for the course-grained representation of the E. coli cytoplasm. B: Weight of each sub-volume for the estimation of the D.

Figure 4 .

 4 Figure 4.19: (A) Example of a QENS spectrum measured on the neutron backscattering spectrometer IN16b at ILL [102] on E. coli cell pellets suspended in D 2 O. The displayed data were recorded at 306 K and q = 0.72 Å 1 . We present the resulting fit of the data assuming three main diffusive contributions accounting for the global and internal motions of proteins, and for the D 2 O dynamics (see section 4.2.2 for further details on the model and the fit procedure). (B-C) HWHM of the Lorentzian components accounting for the global (γ G ) and the local diffusive protein motions (γ L ) at four different temperatures. The dashed lines are the fits of the linewidths assuming the Singwi and Sjölander jump-diffusion model [155]. (D) Apparent self-diffusion coefficient, D G , of an average protein in the E. coli cytoplasm as a function of temperature. A clear non reversible slow-down of D G is visible after the temperature of cell death (T CD ) due to the gelation of the cytoplasm. The inset shows the residence time τ G for the global motions as function of temperature that undergoes an important upturn after T CD and continues to increase during the cooling. (E) Temperature dependence of the diffusion coefficient D L for the local motions of the side-chains. The transition at T CD , where D L starts to increase more steeply with the temperature, is reversible. The inset shows the residence time τ L for the internal dynamics which is nearly constant (⇡ 180ps). (F-G) Geometries of the local motions as derived by the elastic incoherent structure factor (EISF) [52, 55]: p L is the fraction of atoms appearing fixed on the accessible time scale (F), and r L is the radius of the confinement region for this type of fluctuations (G).

Figure 4 . 20 :

 420 Figure 4.20: Comparison of the percentage of H-atoms that appear fixed with those of H-atoms that diffuse in a spherical volume or undergo a 3-jump diffusion.

Figure 4 .

 4 Figure 4.21: (A) Protein diffusion in the CG model cytoplasm simulated by LBMD. Distribution of the traslational diffusion coefficients, and its dependence as a function of the proteins molecular weight. The solid lines represent the fits of the experimental data reported in Ref. [161]. (B) Distribution of local protein concentrations in 17 nm cubic sub-boxes randomly placed in the large simulation system. The vertical lines correspond to the local concentrations of the sub-boxes then back-mapped at the atomistic resolution. (C) Pictorial representation of the CG cytoplasm system (cubic box of side 40 nm) and the schematic strategy of the back-map, and of the temperature scans for the folded and unfolded versions of the atomistic systems. (D) Average protein translational diffusion coefficients, computed in the atomistic systems for the 0.3-5 ns regime and corrected for the effects of periodic boundary conditions, as described in section 4.2.9. The plot shows results obtained with the Amber99SB-disp force field; the results for CHARMM36m are summarized in figure 4.25. (E) Average number of neighboring proteins per protein molecule as a function of its molecular weight. (F) Average number of different classes of protein-protein contacts per protein molecule: hydrophilic-hydrophilic (Phi-Phi), hydrophilic-hydrophobic (Phi-Pho), and hydrophobic-hydrophobic (Pho-Pho). An atom was considered "hydrophobic" if its partial charge was less than 0.2 • e in magnitude. The plots in (E-F) show results obtained from two Amber99SB-disp trajectories, each extended to 1 µs, for the back-mapped all-atom systems with a protein concentration of 288 g/L and simulated at T = 330 K before and after unfolding. The results for CHARMM36m are reported in figure4.25.

Figure 4 . 22 :

 422 Figure 4.22: Proportion of initial protein neighbors of an average protein that remain in contact with the protein in the course of the LBMD trajectory. Two proteins were considered in contact if the shortest distance between their beads was below 7.34 Å.

Figure 4 . 23 :

 423 Figure 4.23: Number of protein clusters present in the course of a sub-volume simulation and the resulting partitioning of proteins into clusters of different sizes for three different protein concentrations in folded (left) and unfolded (right) sub-volumes. The trajectories were obtained with the Amber99SB-disp force field at T = 300 K. Two proteins were considered in contact if the minimum distance between their atom was lower than 3 Å.

Figure 4 .

 4 Figure 4.24: (A) Concentration-dependent diffusion slow-down from simulations. The slowdown is expressed relative to the diffusion coefficient for the least concentrated sub-box (134 g/L). The dashed lines represent fits with a second-order polynomial fraction. (B) Snapshots showing folded-(upper row) and unfolded (lower row) simulation boxes at different protein concentrations. (C) Snapshots of sub-boxes (288 g/L) with a progressively increasing unfolded fraction r u . Folded proteins are blue while unfolded proteins are shown in red. (D) Translational diffusion coefficients in sub-boxes (288 g/L; see panel C) with a varying fraction of unfolded proteins. The insets show the dependence-fitted to a power law-of the apparent unfolded fraction a u on r u . (E-F) The decrease in the translational diffusion coefficient of folded and unfolded proteins inside the partially unfolded sub-box (288 g/L) with increasing r u . The values shown in panels A, E, and F as well as in the inset of panel D are averages across all temperatures, with error bars expressing the standard deviations. The results presented in this figure were obtained with the a99SB-disp force field; analogous plots for CHARMM36m, exhibiting qualitatively the same behavior, can be found in Figure 4.25.

Figure 4 . 25 :

 425 Figure 4.25: Translational diffusion in sub-volumes described using the CHARMM36m force field. (A) Protein translational diffusion coefficients, computed for the 0.3-5 ns regime. (B) Concentration-dependent diffusion slow-down. The slow-down is expressed relative to the diffusion coefficient for the least concentrated sub-box (134 g/L). Averages of the slow-down values for all the temperatures are shown together with their standard deviations. The dashed lines represent fits with a second-order polynomial fraction. (C) Translational diffusion coefficients in sub-boxes (288 g/L) with a varying fraction of unfolded proteins r u ; The inset shows the dependence -fitted with a power law -of the apparent unfolded fraction a u on r u . (D-E) Decrease in the translational diffusion coefficient of folded-and unfolded proteins inside the partially unfolded sub-boxes with increasing r u .

Figure 4 . 26 :

 426 Figure 4.26: Viscosities calculated for the 288 g/L sub-volume at T = 300 K with a varying fraction of unfolded proteins (see Table 4.7).

  Figure 4.26: Viscosities calculated for the 288 g/L sub-volume at T = 300 K with a varying fraction of unfolded proteins (see Table 4.7).

Figure 4 . 27 :

 427 Figure 4.27: Comparison of the transfer function obtained directly from simulations using the Amber99SB-disp force field and a transfer function based on the considerations described in Supplementary Results. Each set of points is fitted with a power law dependence, with the exponents equaling 0.41 and 0.26, respectively.

Figure 4 .

 4 Figure 4.28: (A) Fraction r u of unfolded proteins in the E. coli cytoplasm as a function of temperature. The cell-death temperature of 323.15 K [45], determined from the E. coli growth rate, is indicated by a vertical line. The red solid (Amber99SB-disp) and dashed lines (CHARMM36m) show r QENS u

Figure 4 . 31 : 2

 4312 Figure 4.31: Different components of the MSD for the heating and the cooling phases. ⌦ u 2 ↵ is the total MSD, measured by EINS, taking into account all the types of motions of the average protein in the E. coli cytoplasm. ⌦ u 2 diff ↵ is the MSD due to the roto-tranlational diffusive motions of the entire average protein determined by QENS though the measurement of the apparent diffusion coefficient D G (see eq. (4.26)). ⌦ u 2 int ↵ is the MSD related to the internal sub-diffusive and vibrational motions of the proteins, obtained by the subtraction of ⌦ u 2 ↵ and ⌦ u 2 diff ↵ .

Figure 4 . 32 :

 432 Figure 4.32: Comparison of the internal MSD obtained for the heating with the those calculated for cooling.

Figure 4 .

 4 25 panel C), and this relation also holds for systems with different concentrations (Figure 4.16).

Chapter 5 Differences

 5 between Ca 2+ reach and depleted α-La investigated by MD simulations and NS experiments Based the publication: Differences between calcium rich and depleted alpha-lactalbumin investigated by molecular dynamics simulations and incoherent neutron scattering. Dominik Zeller, Pan Tan, Liang Hong, Daniele Di Bari, Victoria Garcia Sakai, and Judith Peters Physical Review E 101, 032415 -Published 25 March 2020

Figure 5 . 1 :

 51 Figure 5.1: Calcium rich (yellow, bright) and depleted (red, dark) forms of α-La and the variations induced on their structures. α-La from most mammals consists of 123 amino acid residues [169] and its molecular 4 weight is ⇡14.2 kDa.

Figure 5 . 2 :

 52 Figure 5.2: Visualization of the dry and hydrated simulations of α-La dep . (a) Dry environment (0.05h). (b) Hydrated environment (0.4h). The lines show the simulation box size and inside the six chains of α-La dep at 300 K are visualized. The little (blue) molecules indicate water. For the sake of better visibility, we represented the proteins as connected molecules, therefore going beyond the limits of the box, and not as distributed within the same box due to the boundary conditions.

Figure 5 . 3 :

 53 Figure 5.3: Difference of MSD between α-La dep and α-La ca for OSIRIS (a) and SPHERES (b,c). Data are only shown from 100 to 310 K to enhance the differences, since at temperature < 100 K there are no differences within statistical error. The MSD values of the α-La dep are shown as dotted lines and empty symbols, and those for α-La ca as solid lines and filled symbols. The dry samples are blue (lower curves), the 0.4h samples are red (middle curves), and the 0.8h samples are green (upper curves).

Figure 5 . 4 :

 54 Figure 5.4: MSDs from analyzing data from MD simulations using both indirect (GA approximation: lower curves, PK model: middle curves) and direct calculations (upper curves) for the Gaussian resolution function of SPHERES (a), IN13 (b), and OSIRIS (c). On the left side the dry α-La is shown and on the right side the hydrated protein at 0.4h. Each side is ordered in the same way by increasing temperature (280 and 300 K) and α-La dep is next to α-La ca .

Figure 5 . 5 :

 55 Figure 5.5: Distributions of the direct MSDs. Comparison of distributions of the MSDs at t = 30 ps (OSIRIS, highest curve in the peak), 60 ps (IN13, middle curve in the peak), and 1000 ps (SPHERES, lowest curve in the peak) for α-La ca (b) and α-La dep (a). The MSD values were obtained by the average value of the four independent slices of 5 ns for each simulation. ∆ dir in the legend shows the mean value of the distribution as defined in Eq. (5.7).

Figure 5 . 6 :

 56 Figure 5.6: Comparison between the three models obtained by fitting the experimental data (open symbols, dashed lines) and the MD simulations (filled symbols, dotted lines). The circles designate the GA approximation and the squares the PK model.

Figure 6 . 1 :

 61 Figure 6.1: Illustration of the APO and HOLO forms of the MBP. The formation of the MBP+MALT complex lead to an important conformational change, that can be measured by the variation of the open anlge.

Figure 6 . 2 :

 62 Figure 6.2: Dynamic structure factors (left) and vibrational density of states (right) for the MBP+HW, MBP, and HW in the APO and HOLO states.

Figure 6 . 3 :

 63 Figure 6.3: Vibrational density of states (right) for the MBP+HW, MBP, and HW in the APO and HOLO states, calculated by NMA. (left) Estimation of the dynamical structure factor starting as g(E)/E 2 by NMA.

Figure 6 . 4 :κ 1 i 1 σ

 6411 Figure 6.4: (1 st row) Vibrational density of states (vDOS) of selected residues of the MBP.(2 nd row) Order parameters of HW -all the Hydrogen-Bonding Donor/Acceptor, i.e the all the Oxygens, Nitrogens and Sulphurs of the systems (water + protein + maltose) which were less than 4.5 Åaway from water Oxygens were counted.

1

  Note: the equation (3) in the paper of Niessen et al. (http://dx.doi.org/10.1016/j.bpj.2016.12.049) correspond to the projection of the MSPF into an arbitrary direction, i.e.:

  g

  β • r C r β r N r β • r C r β #+ (6.25)An angle of around 160 correspond to the APO form (complexed) of the MBP, meanwhile a value of DOS of approximately 135 , correspond to the HOLO form (uncomplexed) of the MBP. As a reference for this section see(Stockner et. all, 2005).

Table 1 . 1 :

 11 Classification of bacteria according to their optimal growth temperatures.

	Classification	Examples	Optimal Growth Temperature
	Psychrophiles	Psychrobacter articus	T OG < 24 C
	Mesophiles	Escherichia coli (E. coli)	24 C < T OG < 50 C
	Thermophile	Thermus thermophilus	50 C < T OG < 80 C
	Hyperthermophile	Aquifex aeolicus	T OG > 80 C

Table 2 .

 2

	Cold	1 -120	0.09 -10.3	0.13 -1.41	30.8 -3.98
	Thermal	60 -1000	5.17 -86.2	0.99 -4.06	3.98 -0.97
	Hot	1000 -6000	86.2 -517	4.06 -9.95	0.97 -0.40

1: Temperature, energy, and wavelength ranges for neutrons slowed down by different moderators. Moderator Temperature [K] Energy [meV] Speed [km/s] Wavelength [ Å]

Table 2 . 2 :

 22 Coherent, incoherent and absorbtion cross sections (σ coh , σ inc , σ abs ) of the principal elements present in biological systems. The cross section is given for the isotopes of hydrogen and deuterium separately. For the other elements the isotope average is given. All the σ are measured in barns, i.e. 10 24 cm 2 .

	H	D	C	O	N	S	Na	Cl	Al	V
	σ coh 1.758 5.59 5.55	4.23 11.0 1.02 1.66 1.52 1.50 0.02
	σ inc 80.27 2.05 0.001 0.0008 0.5 0.007 1.62 5.3 0.009 5.08
	σ abs 0.333 0.001 0.004 0.002 1.9	0.53 0.53 33.5 0.23 5.08

Table 4 . 1 :

 41 Composition of the average E. coli bacteria (from Neidhardt et al. [104])

	MOLECULES	PERCENTAGE OF DRY WEIGHT
	Macromolecules	96
	Proteins	55
	RNA	20.5
	Lipids	9
	Polysaccharides	5
	Lipopolysaccharides	3.4
	DNA	3.1
	Monomers	3
	Sugars and precursors	2
	Amino acids and precursors	0.5
	Nucleotides and precursors	0.5
	Inorganic ions	1

Table 4 . 2 :

 42 Average parameters resulting from the fit of the Vanadium data with eq. (4.3).

Table 4 . 4 :

 44 Protein composition of the LBMD simulation together with the LBMD coupling coefficient γ for each protein.

	Protein PDB-Structure	MW [kDa] Count γ [10 3 fs 1 ]
	1 TufA	1DG1	43.3	67	1.1920
	2 MetE	SWISS 4zty.1.A	84.5	39	0.6950
	3 CspC	SWISS 3i2z.1.B	7.3	13	3.8150
	4 IcdA	1P8F	91.5	8	0.6478
	5 Pgk	1ZMR	41.1	5	1.1665
	6 FusA	4V9O	77.4	4	0.7346
	7 Eno	1E9I	91.0	3	0.6504
	8 GlyA	1DFO	90.6	3	0.6526
	9 IlvC	1YLR	54.1	3	0.9580
	10 Adk	1AKE	23.6	3	1.4770
	11 Efp	6ENU	20.6	3	1.8653
	12 GapA	1S7C	142.1	2	0.4600
	13 Ppa	2EIP	117.4	2	0.5357
	14 SerC	1BJN	79.3	2	0.7217
	15 DnaK	5NRO	69.0	2	0.8008
	16 CysK	5J43	68.7	2	0.8031
	17 Mdh	2PWZ	64.7	2	0.8400
	18 Tsf	1EFU	60.6	2	0.8815
	19 Tig	1W26	48.2	2	1.0413
	20 SodA	1D5N	45.9	2	1.0779
	21 Upp	2EHJ	45.1	2	1.0928
	22 SodB	1ISC	42.3	2	1.1438
	23 GpmA 1E59	28.4	2	1.5070
	24 Hup	2O97	18.8	2	1.9797
	25 Suc	1JLL	142.0	1	0.4628
	26 PurA	1ADE	94.4	1	0.6326
	27 Asd	1BRM	80.0	1	0.7167
	28 Fba	5VJE	78.0	1	0.7305
	29 PurC	2GQR	54.0	1	0.9591
	30 TpiA	4IOT	53.9	1	0.9597
	31 RpiA	1KS2	45.5	1	1.0860
	32 UspA	SWISS 1jmv.1.A	31.9	1	1.3938
	33 Frr	1EK8	20.6	1	1.8626
	34 PpiB	2RS4	18.2	1	2.0213
	35 SOD1	4BCZ	11.0	10	2.7740

Table 4 . 5 :

 45 Overview of an equilibration protocol used for an initial relaxation of isolated protein geometries. k BB and k SC -force constants of harmonic restraint potentials applied to the heavy atoms of the backbone and the side chains, respectively. ∆t -time step.

	Step Length [ns] ∆t [fs]	Restraints k BB	h	kJ mol nm 2	i	k SC	h	mol nm 2 kJ	i
	1	0.1	1 heavy atoms		1000			500
	2	0.1	1 heavy atoms		500			250
	3	0.1	1 heavy atoms		250			100
	4	0.5	2 heavy atoms		100			50
	5	0.5	2 heavy atoms		50			0
	6	0.5	2	none			0			0

Table 4 . 6 :

 46 Protein concentration and species composition of the five cubic sub-volumes extracted from the LBMD trajectory.

	Conc. [g/L] Protein List
	134.0	CspC (x4), TufA (x3), MetE, Mdh, FusA, SOD1
	197.5	TufA (x4), MetE (x2), CspC, GlyA, CysK, Mdh, SOD1
	279.8	TufA (x2), GapA, GpmA, MetE (x3), CspC, CysK, SodA, IlvC,
		DnaK, Frr, SOD1 (x3)
	287.7	TufA (x3), MetE (x4), IcdA (x2), CspC (x2), PpiB, Mdh, IlvC,
		Efp, GpmA
	297.7	TufA (x8), MetE, CspC (x2), IlvC(x2), SOD1 (x2), Ppa, GlyA,
		Adk, GpmA, RpiA

  .23).

	Preparation of sub-volume geometries
	at different temperatures
	Back-mapped
	configuration
	INITIAL RELAXATION	Folded sub-volume at T i = 300 K
	HEATING PROTOCOL	UNFOLDING PROTOCOL
			T = 300 K
			Time step: 2 fs
		COOLING PROTOCOL	Trajectory length: 10 ns
	Folded sub-volumes	Unfolded sub-volumes
	at T i = 310, 320, …, 380 K	at T i = 300, 310, …, 380 K

Table 4 . 7 :

 47 Protein composition of partially unfolded sub-volumes with a varying fraction r u of unfolded proteins.

	r u = 25 %
	UNFOLDED (5):
	CspC, IlvC, MetE, PpiB, TufA
	FOLDED (11):
	CspC, Efp, GpmA,2x IcdA, Mdh, 3x MetE 2x TufA
	r u = 50 %
	UNFOLDED (8):
	CspC, IcdA, IlvC, 2x MetE, PpiB, 2xTufA
	FOLDED (8):
	CspC, Efp, GpmA, IcdA, Mdh, 2x MetE, TufA
	r u = 75 %
	UNFOLDED (11):
	CspC, 2x IcdA, IlvC, 3x MetE, PpiB, 3x TufA
	FOLDED (5):
	CspC, Efp, GpmA, Mdh, MetE
	4.2.8 Calculation of D t and D r

Table 4 . 8 :

 48 Parameters resulting from the interpolation of solvent viscosities.

	Force field	B w [-] D w [K 1 ] ∆E w [kJ/mol]
	a99SB-disp	20.05	0.0115	24.3
	CHARMM36m 14.75	0.0054	12.8

Table 4 . 9 :

 49 Parameters describing an analytical dependence derived from the computed viscositites for the crowded protein solution.

	Folded
	Force field	B [-] D [K 1 ] ∆E [kJ/mol]
	a99SB-disp	-20.05 32.8272	7157.1
	CHARMM36m	-14.75 12.9697	44982.1
	Unfolded
	Force field	B [-] D [K 1 ] ∆E [kJ/mol]
	a99SB-disp	-20.05 35.2641	36692.3
	CHARMM36m 113942.65 189.9209	251039.8

Table 4 .

 4 10: Different estimations for the fraction of unfolded proteins r u (T ) in the E. coli cytoplasm at increasing temperatures and temperature T 50% needed to unfold 50% of the proteins.

	T 50%

Models

r u (T CD ) r u (T CD + 3K) r u (T CD + 5K)

Table 5 . 1 :

 51 Instrument resolution FWHM vs time. Relation of the FWHM of Gaussian instrument resolution in energy space to the time τ in time space.

	Instrument FWHM (µeV) τ FWHM (ps)
	OSIRIS	24.8	30
	IN13	10.8	60
	SPHERES	0.62	1060

  if we assume that the motion is confined in space and stationary such that: ⌦ r 2 (t + t 0 ) ↵

														t 0	=	⌦	r 2 (t 0 )	↵ t 0	=	⌦ r 2 ↵	(6.12)
	we obtain for the MSD the following expression:
								MSD(t)=2	⇣ ⌦ r 2 ↵	hr(t + t 0 ) • r(t 0 )i t 0	⌘	(6.13)
														3	κ 2	(6.6)
											⌦	u 2 κ	↵	=	1 3	⌦ u 2 ↵	(6.7)
	where	⌦	u 2 ↵	correspond to the Mean-Square atomic Position Fluctuations (MSPF):
		⌦	u 2 ↵	=	⌦	[r(t) l] 2 ↵	⌘	⌦	[r(t) hr(t)i] 2 ↵	=	⌦	r 2 ↵	hri 2 = MSPF	(6.8)
	In this case, the Debye-Waller factor and the MSPF are:
									2W (κ)=	1 3	κ 2 MSPF	(6.9)
							MSPF =	3 2 M h	Z ω m 0	dω	g(ω) ω	• coth	✓	h ω 2 k b T	◆	(6.10)
														t 0	(6.11)

Finally, takinkg into account the atomic Mean-Square Displacements (MSD) defined as:

MSD(t)= ⌦ [r(t + t 0 ) r(t 0 )] 2 ↵

It must also be observed that, in thermophilic bacteria, not only proteins are less affected by high temperatures, but also ribosomes and membranes are more stable[START_REF] Denver | Microbial susceptibility and resistance to chemical and physical agents[END_REF].

Light atoms are generally preferred as moderators since they uptake a high amount of the neutron energy during each collision. Often either H

O or D 2 O is used, where mainly the mass of the H or the D atoms accounts for the moderation efficiency and only about 18 (H 2 O) or 25 (D 2 O) collisions are needed to obtain meV energies[START_REF] Bée | Quasielastic neutron scattering[END_REF].

b α is independent of the neutron energy and is a complex number: the real part can be negative or positive depending of the attractive or repulsive nature of the interaction; the imaginary part represents absorption.

The approximation to the q 2 -term is equivalent to assume that the motion of any atoms in the system can be described by diffusive motion in a harmonic potential that we consider isotropic.

For example, dealing with very light atoms or molecules (e.g. He, H

,D 2 ) or with vibrational motions with characteristic energy comparable or larger than k B T , quantum effects became not negligible.

In MD programs the phase of the simulation used to generate the ensemble is usually named: equilibration.

For biological MD simulations ∆t is usually of the order of few femtoseconds (10 15 s) therefore, in order to obtain a trajectory of few nanoseconds (10 9 s), one has to perform at least a million of integration steps.

A force field, in the context of a computer simulation, refers to the functional forms used to describe the intra-molecular and inter-molecular potential energy of a collection of atoms, and the corresponding parameters that will determine the energy of a given configuration. Thus it is a special case of interatomic potentials and it must not be confused with force field in classical physics.

Indeed, the number of covalent bound is proportional to the number of atoms.

In the article: Scalable Molecular Dynamics with NAMD (Journal of Computational Chemistry, 2005), Phillips and his colleagues note that this method, implemented in the NAMD software, is symplectic and time reversible[START_REF] Phillips | Scalable Molecular Dynamics with NAMD[END_REF].

[START_REF] Musto | Genomic GC level, optimal growth temperature, and genome size in prokaryotes[END_REF] When good energy conservation is needed for NVE ensemble simulations, it is recommended to choose 2 fs,2fs,and4fs as the inner, middle, and outer time steps if rigid bonds to hydrogen atoms are used; or 1 fs,1fs, and 3 fs if bonds to hydrogen are flexible. More aggressive time steps may be used, instead, for NVT or NPT ensemble simulations -i.e. 2 fs,2fs, and 6 fs with rigid bonds and 1 fs,2fs, and 4 fs without[START_REF] Phillips | Scalable Molecular Dynamics with NAMD[END_REF].

Usually, in this phase of the simulation, periodic boundary conditions are particularly important to avoid surface effects.

This might be done also in alternation with some phases of cooling.

Clearly, between the velocities of the atoms and the temperature of the system there is a strong connection. According to the classical equipartition theorem, each normal mode has[START_REF] Dill | Driving forces in the origins of life[END_REF] 2 k B T energy, on average, at thermal equilibrium. Thus hEi =[START_REF] Dill | Driving forces in the origins of life[END_REF] 2 ∑ i m i v 2 i = 1 2 Nk B T . One of the most common ways to set the initial velocities is by generating numbers pseudorandomly and choose them so that the total kinetic energy of the system corresponds to the expected value at the target temperature T .

Clearly, an exception is made for the amino acid proline, whose the side chain is formed by only one H atom. Therefore it is represented in the OPEP model by all its atomic constituents.

Further details on the connection between the apparent and real unfolded fraction are reported in section 4.4.2.

Acknowledgments D.Z. was supported by a Ph.D. scholarship cofunded by the Communauté Université Grenoble Alpes, the STFC Rutherford Appleton Laboratory, and the Institut Laue Langevin. D.D.B. acknowledges support via a scholarship from the da Vinci program of the French-Italian University. J.P. and V.G.S. gratefully acknowledge the support by M. Johnson to partly acquire the Ph.D. grant. The authors gratefully acknowledge the financial support provided by JCNS to perform the neutron scattering measurements at the Heinz Maier-Leibnitz Zentrum (MLZ), Garching, Germany. Experiments at the ISIS Neutron and Muon Source were supported by a beam time allocation from the Science and Technology Facilities Council. We thank the ILL for beam time allocation. D.Z. thankfully appreciates the assistance of M. Zamponi during the experiment at the MLZ and the MDANSE software support of E. Pellegrini and R. Perenon.

Γ = 86. For sake of comparison, we recall that Dill obtained the following values ∆H † ' 27 kJ/mol and Γ = 51 but starting with a very different form of the stability curve. , and the temperature dependent diffusion constant measured by QENS. The fit is performed using a constant intrinsic growth-rate parameter g 0 (blue line), a temperature dependent g 0 (T ) ∝ D QENS (orange), and assuming a negligible activation barrier ∆H † (green).

We then introduce into the fit the temperature dependence of the exponential prefactor in terms of a reaction-diffusion model by assuming g 0 ∝ D QENS G (see Fig. 4. [START_REF] Tsuchido | Destruction of the outer membrane permeability barrier of Escherichia coli by heat treatment[END_REF]). The fit, see Fig. 4.29, is comparable to the case where g 0 is constant, but by including the dependence on diffusivity we recover smaller activation barrier and number of essential proteins (∆H † ' 35 kJ/mol and Γ = 81). A final numerical test is done by assuming the growth rate of the bacterium completely rate-limited by diffusivity (∆H † = 0). In this case, the E. coli growth rate curve cannot be fitted at temperatures below the cell-death.

The obtained results confirm that even without assuming a proteome catastrophe, the growth-rate of E. coli can be reproduced very well by a more smooth temperature progress of unfolding of the proteome; and that including the diffusion contribution to reactivity may help tuning the essential parameters of the model, the activation barrier and the estimate of essential proteins for the bacterium growth.

Appendix

EINS Results

The temperature dependence of the total MSD measured with EINS is similar to the one observed for the apparent diffusion coefficient D G measured by QENS and the parameters T 0 and ∆T resulting from the fit with eq. The total MSD represent an average over various motions in the time limit going to infinity. However, as described in sections 2.2.3 and 2.2.2, typical biosystems contain several kind of motions and following [START_REF] Bée | Quasielastic neutron scattering[END_REF][START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF], we can subdivide the movements into a sum of three contributions: local vibrations, diffusive motions of molecular sub-units and global diffusion of the entire protein, as they correspond to distinct time scales. Consequently, we can rewrite the total MSD, ⌦ u 2 ↵ , as the sum of three different MSD:

On the other hand, it is possible to show that ⌦ u 2 diff ↵ can be derived from the global diffusion coefficient [START_REF] Hennig | Dynamics of highly concentrated protein solutions around the denaturing transition[END_REF]:

where t res is the resolution time determined from the energy resolution ∆E FWHM of the instrument [START_REF] Zorn | On the evaluation of neutron scattering elastic scan data[END_REF].

We therefore combined EINS and QENS to measure

However, the observed dynamical transitions is not reversible-after the transition, D G (T ) remains equal to D 2 (T ), even at low temperatures. The same dynamical behavior was previously found for proteins in solutions and, due to the similarity of the systems, this suggests that the underlying process is the same. In particular, in the first state (state 1), the E. coli cytoplasm is a liquid where the majority of proteins are folded. After the transition (state 2), the cytoplasm is in a gel state, and the proteins are mainly unfolded. For simplicity, we will assume that in states 1 and 2, all the proteins are fully folded and fully unfolded, respectively. Therefore, the diffusion coefficient D 1 will describe the global diffusion of the folded average protein in the liquid state, termed D ( f ) , whereas D 2 will be related to the global diffusion of the fully unfolded average protein in the gel state, termed D (u) . Thus, for D G (T ) we have:

In this picture, a QENS u (T ) can be interpreted not only as a measure of the diffusion slow-down during the unfolding, but it also represents the apparent fraction of unfolded proteins that weights D ( f ) and D (u) and that is necessary to reproduce the values of D G for systems in a partially unfolded state. In particular, our simulations suggested that the relation between the apparent fraction of unfolded proteins determined dynamically, a MD u , and the real fraction of unfolded proteins, r MD u , can be described by a power law:

Therefore, to estimate the real unfolded fraction as a function of temperature starting from the dynamical QENS results, we can use the inverse relation:

In conclusion, comparing QENS experiments and MD simulation, first we verified that the simulations can reproduce the linearity of D 2D and4.25A). Then, since protein unfolding is computationally demanding for all-atom MD simulations, instead of calculating extremely long trajectories to allow all proteins to unfold, we forced the unfolding (as described in Subsection "Preparation of partially unfolded sub-volumes"), and we performed a series of relatively short simulations (⇠100 ns), sufficiently long to determine the diffusion coefficient in the timescale explored by the QENS experiments, but brief Supporting Information [START_REF] Gerald R Kneller | nMOLDYN: a program package for a neutron scattering oriented analysis of molecular dynamics simulations[END_REF]. Both methods use the 20 ns trajectories.