The bending of hemagglutin

(1)

The bending of hemagglutin

A bachelor research about the bending of the hinge region of hemagglutinin and enhanced sampling in molecular dynamics via

hamiltonian replica exchange with solute tempering.

by Cees de Wit

Supervised by:

S. Boonstra, MSc Prof. dr. ir. E. van Giessen

Prof. dr. ir. P. R. Onck

Rijksuniversiteit Groningen

July 12, 2016

(2)

1 Abstract

When a influenza virus tries to enter a cell, it has to fuse its viral membrane with the cell membrane. This fusion is mediated by a protein called hemagglutinin. hemagglutinin(HA) facilitates the fusion by bending into a hairpin configuration. In this study the bending region, called the hinge region, was simulated to elucidate the molecular dynamics(MD) and structures of this hinge region. To enhance the MD simulations, Hamiltonian replica exchange was used. Two implementations of Hamiltonian replica exchange where tested, a implementation using the free energy perturbation of GROMACS and a implementation using the plugin PLUMED. The methods where tested on an alanine dipeptide and a 3K(I) molecule. The PLUMED method was easiest to implement and had the fastest simulations, thus was used for the HA simulation.

The hinge region was simulated from an initial straight, helical configuration. The histidine residue of the peptide was not protonated. After about 700 ns of simulation the peptide was in a folded configuration, in contrast to simulations done by Kalani et al..

The configuration switched between a helical folded structure and a non helical folded structure. In the folded configuration, the location of the histidine residue was on the outside of the folded structure, the location of the non polar residues was on the inside and the acid, polar and non polar residues on the outside of the structure. I propose that the bending of the hinge region is probably mediated by hydrophobic interactions and not by the protonation of the histidine residue.

2 Introduction

It is winter, you wake up and feel horrible. You have a headache, are nauseous and are burning up. You have the flu. It is a feeling that everybody has experienced and it is the result of a tiny virus, called the influenza virus. In my bachelor research I modulated a protein of the influenza virus, called hemagglutinin. To do this I used a computational method, called molecular dynamics. In the next sections some background information about the influenza virus and hemagglutinin is given. Then, a method to enhance the molecular dynamics simulations, called replica exchange, is explained. Two variants of the replica exchange were tested and both will be explained. Thirdly, the details of the simulations that were done and their results are discussed. The report ends with some conclusions about the simulations. I hope you enjoy reading this report.

(4)

3 Physical background

3.1 Influenza

The virus that gives us the flu is called the influenza virus. Influenza is an envelope virus, which means that its RNA is encapsulated by a lipid bilayer. To infect a cell, the virus has to release its RNA into that cell. Then, the RNA translated by the host cell and the virus is replicated.

Figure 1: Endocytosis of the influenza virion [17]

An infection starts with the entering of a virus into the cell via endocytosis. The process of endocytosis is shown in Figure 1. After endocytosis the virus is in the cell, but it is still separated from the cytoplasm by a cell membrane. To release its viral RNA into the cell, the viral membrane has to fuse with this cell membrane. The fusion between the membrances is thermodynamically favorable, but first a kinetic barrier must be overcome. This kinetic barrier is created by repulsive hydration forces between the membranes and its repulsion increases steeply when the two membranes come within 2 nm from each other.[12] For several viruses, e.g. the influenza virus, this kinetic barrier is overcome by a protein called hemagglutinin.

3.2 Hemagglutinin

Hemagglutinin (HA) is a protein that is located on the outside of the viral membrane.

A large fraction of the parts sticking out of the virus in Figure 1 are HA proteins. HA mediates the fusion between the viral and cell membrane. This is done in several steps, which are depicted in Figure 2.

(5)

Figure 2: Schematic of membrane fusion mediated by HA. The upper membrane is the cell membrane and the lower is the viral membrane. In step 1 and 2 the whole protein is shown. In step 3, 4 and 5 only HA₂ is depicted. [12]

The hemagglutinin consists of two main parts, called HA₁ and HA₂. In Figure 2, HA₁ is red and HA₂ is blue. The main function of HA₁ is to recognize a proper host cell and

”dock” to it . The main fuction of HA2 is facilitating the fusion.[13] In the first step a receptor binds to HA₁. In the second step, the HA₁ heads separate as a result of a reduced pH in the endosome. In step 3, an extended intermediate is formed. The HA₂ part is stretched out between the two membranes and forms a kind of bridge between them. Then, the actual fusion begins in step 4. The extended HA₂ folds together, like a hinge, which pulls the two membranes to each other. After the bending, the fusion pore is formed and stabilized, this is shown in step 5. Now, a pore is created and the viral RNA can be sent into the cytoplasm.

In the next section, we will go into more detail of the folding in step 3 and 4, which was the focus of the simulation done.

3.3 The hinge peptide

It remains elusive what the dynamics of the configurational rearrangements in step 3 and 4 of Figure 2 are, as well as how a lower pH exactly influences the folding. In a research done by Kalani et al. the effect of a lower pH on a part of the HA₂ peptide was investigated via molecular dynamics(MD) simulations. This research is discussed in the next part.

In Figure 3, a fusion scheme of HA is shown again. HA₂ consists of three identical monomers. In B1, B2 and B3 one such monomer is depicted. The monomer is shown in its extended configuration, connecting both membranes, in B2. In B3 the folding occurs.

The HA₂ monomer bends at a certain region called the hinge region. It was this region that Kalani et al. and I simulated.

(6)

Figure 3: Schematic of membrane fusion mediated by HA. In A, the whole protein is shown, the blue part on the top is HA₁, the red and purple part is HA₂. [13]

In the Kalani paper it was hypothesized, that the conformational changes due to the low pH in the endosome were mediated by histidine residues. This is because histidines are the only residues in the protein that have a pKa matching the pH of the late endosomes(pK_a of 6.)[14]. This means that only the histidine is protonated during the pH drop.

To test how the protonation of a histidine influences the bending of the hinge region, Kalani et al. did several MD simulations. Simulations were done using residues 95 to 120 of the HA. This peptide consisted of the hinge region, 6 residues large, flanked by two arms of a length of 10 residues. This peptide will now be referred to as the ”hinge peptide”.

In their study, they used several hemagglutinin types. For each type they did simulation in an acid environment and a neutral environment. To simulate an acid environment, the histidine in the hinge region was protonated. The simulations were started from a straight configuration, mimicking step B2 in Figure 3, or a bent configuration, mimicking step B3. Some simulations were done using an implicit solvation system and some using an explicit solvation system.

The results they got from their simulations were as follows.

- For unmutated systems, all simulations running at neutral acidity ended with a straight configuration, no matter if the starting configuration was bent or straight - For the acid simulations, all unmutated systems ended with a bent configuration,

no matter the starting configuration.

- For systems with a mutation of the histidine to an alanine, all end configurations were straight.

From the results of the simulations the authors propose, that the protonation of the Histidine in the hinge region may be of importance to the bending of HA.

(7)

These results are pretty clear, but the simulations were done using a CHARMM27 force field.[7][6]. This force field is known to result in conformations with a higher amount of helicity then in experiments.[3] This may have lead to structures with too much helicilty for the neutral simulations, which could have lead to the straight end configurations.

In an updated version of CHARMM, called CHARMM36 [3], this tendency to helical structures was lessened. In the paper also regular MD was used. It may be that the simulations got stuck in a configuration, which did not have the lowest energy. An enhanced sampling method, would maybe find other configurations. In my simulation the updated CHARMM36 force field and enhanced sampling was used, to study if the hinge peptide is in a straight or bent configuration at neutral acidity. In the next section MD and an enhanced sampling method is explained.

4 Molecular dynamics

4.1 Introduction into molecular dynamics

Molecular dynamics (MD) is a computational method to simulate molecules. The program that is used for the simulations, is called GROMACS[23]. MD is often used to simulate biological assemblies, like membranes and proteins. It allows to see how proteins behave in length and time scales that are impossible to see in e.g. microscopes.

In a simulation, the forces on atoms are calculated and used to determine the next configuration of a molecule. Molecules can have many different configurations with many different energies. A map of all the possible conformations and their energies is called the energy landscape. A example of a energy landscape is shown in Figure 4. One of the challenges of MD is to find the lowest energy configuration in the energy landscape.

The energy landscapes of proteins often have several local energy minima, separated by energy barriers. When running MD simulations, there is a change that a simulation gets stuck in such a minimum. A method to overcome this, is to make use of replica exchange, which is explained in the next section.

Figure 4: Energy landscape of a protein [19]

(8)

4.2 Theory of normal replica exchange

With replica exchange molecular dynamics (REMD) several replicas of the simulation are made. The replicas usually differ from each other in their temperature, but can also differ in, for example, their Hamiltonian. For clarity, temperature exchange will be used in this explanation. The lowest energy replicas correspond to the original simulation, while higher replicas have a higher temperature and thus a higher energy. The replicas with a higher temperature may have enough energy to get over the energy barriers and sample more of the energy landscape. Conformations between different replicas are exchanged, such that the lowest replica can get the conformations of the higher energy replicas.

The replicas are exchanged via a Monte Carlo algorithm[2]. For temperature REMD the probability that a replica is exchanged is given by

P (1 ↔ 2) = min

1, exp

1

k_BT₁ − 1 k_BT₂

(U₁− U₂)

, (1)

where T1 and T2 are the temperatures and U1 and U2 are the instantaneous potential energies of replicas 1 and 2 respectively. Here, replica 1 has a lower temperature then replica 2.[23]

A problem with REMD is that the number of replicas increases with the square root of the degrees of freedom. For large systems this becomes computationally very expensive, because more and more replicas need to be simulated.[20] Because of this, it is desirable to decrease the number of used replicas. With temperature REMD the dynamics of the whole system (protein and solvent) is enhanced, while in general we are only interested in enhancing the protein dynamics. If only the protein dynamics is enhanced, the energy difference between the replicas is smaller and fewer replicas are needed. A way to do this is called replica exchange with solute tempering(REST)[20] and is described in the next section.

4.3 Theory of replica exchange with solute tempering

As mentioned before, with REST only the dynamics of the protein is enhanced. In the MD program GROMACS, this can be done via a rescaling of the force field parameters.

An implementation of a variant of REST that was proposed by Terakawa et al was used in this study.[20].

The idea is to give the protein a different ”effective temperature” for every replica while the water temperature is the same for all the replicas. This can be done via scaling of the Hamiltonian. The Boltzmann distribution contains the temperature and the Hamiltonian in the form βH in which β = 1/kBT . So thermodynamically speaking, reducing the Hamiltonian corresponds to increasing the temperature.[20][8] When replicas with different Hamiltonians are made the same Monte-Carlo algorithm is used as for regular temperature REMD only with a different exchange probability, namely

P (1 ↔ 2) = min

1, exp

1 k_BT

(U₁(x₂) − U₂(x₁)) + (U₁(x₁) − U₂(x₂))

(2) where T is the temperature of the replicas and U_i(x_j) are the potential energies of replicas i with the coordinates from replica j.

(9)

To implement REST, different replicas have to be made. First the potential energy function must be divided in a solvent part, a protein part and a protein-solvent part.

This can be done by dividing the potential energy function three terms as

V = V_pp+ V_ps+ V_ss, (3)

in which V_pp is the protein-protein, V_ps the protein-solvent and V_ss the solvent-solvent interaction potential.

Using this expression two end potentials can be defined. One end potential corresponds to the potential of the system at the temperature of interest and one to the potential of the replica with the highest ”effective temperature”. The potentials are as follows,

V_L= β_L β V_pp+

s β_L

β V_ps+ V_ss (4)

V_H = βH

β V_pp+ s

βH

β V_ps+ V_ss (5)

where β is the inverse of the system temperature and β_L and β_H are the inverse of the lowest and highest effective temperature respectively. Normally β_L = β. The factor β_H/β reduces the protein-protein interactions, thus increasing the solute dynamics. Using his potential the intra protein interactions are scaled by β_H/β, the protein-solvent interaction bypβ_H/β and the solvent interactions are not scaled. In an earlier version of the method, called REST1, the arithmetic average of the scaling for the protein-solvent was taken, instead of the square root. The square root scaling of the protein-solvent allows for quick implementation. This implementation was tested by Berne et al.[24] and is called REST2.

REST2 has better performances than REST1.

These other potentials can now be defined as a linear combination of the two end potentials. The potentials are made using a variable called λ_i, which has a value between 0 and 1. The potential for replica i is defined as follows,

V_i = (1 − λ_i)V_L+ λ_iV_H

= β_L(1 − λ_i) + β_Hλ_i

β V_pp

+

s β_L

β (1 − λ_i) + s

β_H

β (1 − λ_i)

V_ps+ V_ss

(6)

The potentials for λ = 0 and λ = 1 correspond to VL and VH respectively. Other λ⁰s result in a potential that is a linear combination of V_L and V_H. These different potentials will define the different replicas. In the next section will be explained how this method is implemented in GROMACS.

(10)

4.4 REST in GROMACS

GROMACS supports Hamiltonian replica exchange (HREX) by making use of its free energy functionality. In the free energy function of GROMACS, a perturbation can be made between two states, e.g. between ethane (state A) and methane (state B). Then, for example the free energy difference between the two states can be calculated. These perturbations make use of a linear λ function, much like the one in equation 6, between the states. For example, the interpolation of a bonded potential is

V = ((1 − λ)k^A+ λk^B)f, (7)

in which k^A and k^B are the force constants of state A and B and f is the rest of the potential.

A state A can be made in which the potentials correspond to VL and a state B can be made in which the potentials corresponds to V_H. Using these as the different states for the free energy functionality, the HREX of GROMACS can be used to implement REST.

How these states are created is discussed next.

The potential of the system generally consists of the Coulomb interaction, the Lennard- Jones(L-J) potential and bonded potentials. These are often defined as follows,

V = V_c+ V_LJ+ V_bonded

= 1

4π₀ X

i<j

q_iq_j

r_ij +X

i<j

_ij σ_ij r_ij

12

− σ_ij r_ij

6

+X1

2k^b_ij(r_ij − b_ij)²+ ... .

(8)

In this equation, q_i is the partial charge of atom i, ₀ is the dielectric constant, r_ij is the distance between atoms i and j, σ_ij is the average van der Waals radii of the atoms i and j, kij is the spring constant, bij the bond length in equilibrium and ij =√

ij, where i

is the depth of the Lennard-Jones potential well.

Based on the standard potential V, V_H can be constructed by rescaling the parameters in the parameter files of GROMACS. First the scaling factor is defined, which is γ = βH/β.

Then the different potentials are scaled as follows:

- For the Coulomb potential the partial charges of the protein are scaled as q_H,i =

√γqi. This scaling is with a square root, since the charges are multiplied in the coulomb potential. leading to a scaling of γ.

- For the Lennard-Jones term the parameters are scaled as _H,i = γ_i if i is in the protein.

- All the spring constants in the bonded term that belong to the protein are scaled with γ, for example k_ij^b = γk_ij^b.

If all these rescalings are implemented simultaneously, the rescaling of equation (5) is realized.

(11)

4.5 Hamiltonian replica exchange with PLUMED

Because the simulations using the REST implementation were relatively slow, a different method was tried. In this method, a plugin for GROMACS called PLUMED [21] is used.

PLUMED supports an implementation of HREX made by Bussi [8]. This implementation does not use the slow Lambda perturbation of GROMACS, which makes the simulations faster. In the next section the theory of this implementation is explained. Then, the simulations to test PLUMED and their results will be discussed.

What is PLUMED?

PLUMED is a plugin for MD simulation software, like GROMACS, that can perform calculations while a simulation is running and can integrate those with that simulation. In Figure 5, the work flow of the program is shown. The first step occurs before the simulation runs. PLUMED is set up from the details of the simulation and the PLUMED input file. From this step, PLUMED knows what to do during the simulation. In the second step, the simulation really starts. First, the forces are calculated by the MD program.

Then the atom positions are sent to PLUMED and PLUMED does the calculations that were specified in the setup. When the calculations are finished, PLUMED sends the new configuration back to the program. The results are integrated and the cycle begins again with the calculation of the MD forces. This process is repeated until the simulation is finished.

Figure 5: PLUMED workflow [22]

(12)

REST in PLUMED

The implementation of Bussi makes use of a REST2 method, just like the implementation suggested by Terakawa. The main difference between the implementations is thus not how the Hamiltonian is modified, but how it is implemented in GROMACS. Instead of using the Lambda functionality of GROMACS, the HREX functionality of PLUMED is used.

The different replicas are created by editing the topology files. For each replica the force parameters are scaled corresponding to the desired modification of the Hamiltonian.

The parameters are scaled via a scaling factor η. η = 1 corresponds to the replica at the temperature of interest. The other ηs depend on the desired effective temperature. For example, the start temperature 300K and the desired effective temperature 350K give η = 300/350 = 0.8571. Using this η, the force -field terms of the protein are scaled in the following way:

- The charge of the atoms in the protein is scaled with √ η.

- The Lennard-Jones term is scaled by η.

- The proper dihedrals for which the first and the fourth atom are of the protein are scaled by η.

- The atom pairtypes and c-map are scaled by η.

This scaling is so that the interactions in the protein are at an ”effective” temperature of T /η, the interactions between the protein and the solvent at T /√

η and the interactions only in the solvent at a temperature T . The ”effective temperature” does not represent the real temperature, but indicates the increase of the protein dynamics.

With the different replicas created, PLUMED can be used to do the HREX. The HREX is performed in the calculations step depicted in Figure 5. The work flow of this calculation is shown in Figure 6.

(13)

Figure 6: HREX workflow[8]

At the first step an unconditional replica exchange is performed. After this exchange the total energy of the systems is calculated. This energy is stored for a later step.

In the next step, again an unconditional replica exchange is performed, so that the original(unswapped) states are restored. Then the ”real” replica exchange is performed.

For this the stored total energy is used to calculate the acceptance ratio of the exchange used for a Monte Carlo exchange. The exchange probability is calculated using equation (2). The results are sent back to GROMACS were it is implemented in the simulation.

After several simulation steps are done, a new replica exchange is tried and the process repeats itself.

PLUMED has a simple script called, the partial tempering script, that can be used to implement REST. With a small modification, needed to modify the CMAP[16], the script changes all the necessary force field parameters. This script makes the implementation of REST2 via PLUMED quick and easy. In the appendix is explained how this script is applied and how the HREX is performed in GROMACS.

5 Simulation protocol and methods

In this section, the preparation and details of the simulations are discussed. Also the methods to analyze the simulations are explained. In the first section, the settings that are almost the same for each simulation is explained. Then the other details are discussed per simulation.

The two implementations of REST, discussed in the previous sections, will be called REST and PLUMED, even though both methods make use of the REST2 scaling. To see if implementation of REST and PLUMED work in GROMACS, some test simulations were done. For the simulations two peptides were used.

The first peptide was an alanine dipeptide. This peptide was chosen, because it is

(14)

small peptide and a small peptide allows for quick sampling and testing. Another reason to chose alanine dipeptide, is that is tested alanine dipeptide in the Terakawa paper, so the paper could be used as reference material. For the alanine dipeptide, first a normal REMD run was done to serve as reference. Then a REST and PLUMED simulation was done.

The second peptide used, was the 3K(I) peptide. This is a peptide consisting of 13 Alanines and 3 Lysines. 3K(I) can be folded in a specifically stable helical formation[15].

This peptide was chosen because my daily supervisor, Sander Boonstra, researched 3K(I) and did a normal temperature replica exchange simulation using this peptide. The results of that simulation could be used as reference. The REMD simulation could also be used to compare the speed of REST and PLUMED to normal REMD for larger peptides.

Lastly a simulation on the hinge peptide was done using PLUMED. The details of this simulation are described in the fourth section.

5.1 Preparation of the systems

The force field used for the simulations was CHARMM36. The systems were solvated using the TIP3P watermodel. The system was simulated using dodecahedron boundary conditions and the particle mesh Ewald method[11] was used. The systems were equilibrated by first a steep energy minimization on the unsolvated system, then by another steep energy minimization on the solvated system, followed by a particle restraint energy minimization. The equilibration were finished by a constant volume(NVT) and then a constant pressure(NPT) simulation. The NVT and NPT runs were done at a constant temperature, corresponding to the temperatures of the lowest replica. After the NPT simulation the box was scaled, so that the end volume of the box matched the average volume of the NPT run. This is necessary, since the box size fluctuated during the NPT simulation. Because of this, there is a change that the simulation ends in a very small or very large box size. This would result in a high or low pressure. The average volume is used to get the average pressure, a pressure around 1 bar, for the final simulations. This process will be referred to as the energy minimalization .

5.2 Alanine dipeptide

The alanine dipeptide was made using the atom builder PyMol[18]. To an Alanine molecule, an Acetyl group and a Methylamide group was added. The names in the topology were edited so that they matched the names specified in the CHARMM forcefield files. Using the Alinine Dipeptide, three simualtions were done. First the REMD simulation is described.

REMD simulation

For the REMD simulation 32 replicas were used. The temperatures of the replicas ranged from 300K to 460.16K. To create the different replicas, the energy minimalization was done until and including the position restraint run. The box size was so that the simulation included 1621 water molecules. Then the different replicas were created by doing the NVT and NPT runs at different temperatures corresponding to the different replica

(15)

temperatures. Again the boxes were scaled to make sure the pressure was correct, and replicas were ready.

REST simulation

For the REST simulation only 3 replicas were needed. The three replicas were made from the REMD replica at 300K. Then, for each replica the topologies were edited so that the scaling described in paragraph 3.4 was achieved. The replicas were created using the variable λ as in equation (6). The values of λ used were 0, 0.25 and 0.5. Here the highest replica had an effective temperature around 450K.

PLUMED simulation

The PLUMED implementation only works if the vdw and rcoulomb and rlist are the same size. The value of vdw, which was 1.4, was used, because this would most probable have the smallest implications. It means though, that the simulations were slower than if a smaller cutoff length was chosen. To create the replicas again the REMD replica at 300K was used. The different replicas were created using the partial tempering script with a scaling of 1, 0.707 and 0.5. The highest replica had an effective temperature of 600K.

5.3 3K(I)

For the 3K(I) peptide only a REST and PLUMED simulation were done. I used a configuration that I got from Sander, my supervisor. In this configuration the peptide was straight and had no secondary structure.

REST simulation

The energy minimalization was used to prepare the system. The box had 2392 water molecules and the temperature of the lowest replica 275 K For the simulation 5 replicas were used and they were created the same way as the replica of the alanine dipeptide.

The Lambdas used were 0.01, 0.12, 0.23, 0.34, 0.45. The highest effective temperature was around 427K.

PLUMED simulation

For the PLUMED simulation the lowest replica of the REST simulation was used. Then, the partial tempering script with CMAP correction was used with the scaling factors 0.644028, 0.703273, 0.767968, 0.838614, 0.915757, 1, to create the different replicas. Again the highest effective temperature was around 427K.

5.4 Hinge peptide

For the HA simulation a protein structure was used from the Protein database(REF).

The structure ID is 1HGF. From the structure all molecules were deleted except for one monomer. From this monomer only residue number 95 up to 120 were simulated. This is the hinge region with two arms as described in section 3.3. The energy minimalization was done on the peptide at 275 K. The box size was solvated with 5838 water molecules.

Then again, the partial temperingscript with CMAP was used to create the different

(16)

replicas. The scaling factors used were 0.6, 0.64, 0.68, 0.72, 0.76, 0.80, 0.85, 0.95, 1. This simulation ran on the supercomputer Cartesius, on 128 nodes for 700 ns.

5.5 Analysis methods

In this section a few analyze methods, that were used, are explained.

To see how well the REST and PLUMED simulations samples the conformational space, Ramachandran plots were compared. In a Ramachandran plot, two dihedral angels in the backbone of a protein are plotted against each other. The two angels plotted are called φ and ψ and are depicted in Figure 7a. The combination of these angles shows in which form the peptide is, e.g. α helical or in a β sheet. In Figure 7b, the regions corresponding to α helices and β sheets are indicated.

(a) The phi and psi dihedral angles of the alanine dipeptide

(b) Different regions in the Ramachan- dran plot.

Figure 7: Properties of the Ramachandran plot

For analysis of the HA simulation two methods are used. The first method is an angle calculation, as is done in the Kalani paper. The angle was defined as the angle between the two arms of the hinge and is shown in Figure 8. The beginning of the arms was defined by the Cα atoms of the second last residues of the peptide, residue 96 and 119.

As end points, the backbone nitrogen(N) of residue 100 and 115 were used.

(17)

Figure 8: Angle used for the HA calculations. The histidine 106 and aspartate 112 residues are shown explicitly.

The second method was to calculate the root mean square distance, dRM S The dRM S

is used to see how closely the result of the HA simulation resembles a post fusion crystal structure. To calculate the d_{RM S}, the distances d⁰_ij between residues that are in contact in the crystal structure are measured. These contacts, called native contacts, are defined as all N_ij pairs that lie within 9 ˚A of each other in the crystal structure. Pairs that are within 9 residues form each other in the sequence are excluded. For each conformation x in the trajectory, the root mean square deviation of distances dij(x) relative to the native contact distance, is calculated as follows[4]

d_{RM S} =

1 N_ij

X

i,j∈native

(d_ij(x) − d⁰_ij)²

1/2

. (9)

6 Results and Discussion

In this the results of the simulations will be shown and discussed. Fisrt the speed of the simulataions is discussed, which shows why plumed was used. The results for the alanine dipeptide and 3KI simulations are discussed and finally, the results of the HA simulation are discussed.

Speed of the simulations

In table 1, the speed of the alanine dipeptide and 3K(I) simulations are shown. The third and fourth columns show the number of replicas used for the simulations and on how many computer nodes the simulations ran. The fifth and sixth column show the speed of the simulation and the speed of the simulation per node. On remark about the table, the PLUMED 3K(I) simulations is not the one described in section 5, but a run with similar properties, only with 5 replicas. This allows for a better comparison with the REST simulation.

In the table it can be seen that the REST simulations are slower than the PLUMED simulations, especially for the 3KI(I) peptide. The reason that the REST simulation are slow is probably that the free energy perturbation of GROMACS is rather slow. If the peptide is large, more of those free energy perturbation have to be done. How larger the

(18)

Molecule Method # Replicas # Nodes Speed Speed/node

Alanine dipeptide REM 32 32 9 ns/day 9 ns/day/node

Alanine dipeptide REST 3 30 30 ns/day 3 ns/day/node

Alanine dipeptide PLUMED 3 30 44 ns/day 4.4 ns/day/node

3KI REM 32 128 30 ns/day 7.5 ns/day/node

3KI REST 5 30 9 ns/day 1.5 ns/day/node

3KI PLUMED 5 30 18.3 ns/day 3 ns/day/node

Table 1: Speed of the simulations

system, how larger the speed difference between PLUMED and REST becomes. Since PLUMED is faster, it was used for the 3K(I) and HA simulations.

6.1 Alanine dipeptide

From each alanine dipeptide simulation a Ramachandran plot was made. Figure 9 show the plots of the REMD, the REST and PLUMED simulations. Of each simulation there is a plot of the highest and lowest temperature replicas.

From the plots it can be seen that in the high temperature replicas, Figure 9d, 9e and 9f, more conformations are sampled than in the low temperature replicas, Figures 9a, 9b and 9c.

The low temperature plots are almost the same for each simulation. For the high temperature plots the REMD and REST simulation plots look like each other. The difference between those could be addressed to more sampling for the REST simulation.

The PLUMED simulation samples a larger space than the other two. This is because the PLUMED replica has a higher temperature than the REMD and REST replicas.

From the likeliness of the REMD and PLUMED plots to the REMD plot, it can be conclude the methods sample the conformational space well.

(19)

-150 -100 -50 0 50 100 150 -150

-100 -50 0 50 100 150

Psi

Phi Ramachandran Plot

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(a) REMD simulation 300K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(b) Ala REST simulation 300K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(c) Ala REST simulation 300K PLUMED

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(d) REMD simulation 460.1K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(e) Ala REST simulation at 460K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(f) Ala PLUMED simulation at 600K

Figure 9: Ramachandran plots of the alanine dipeptide simulations

The Ramachandran plots do not have an entirely correct shape. even though the simulations give good results compared to each other. In Figure 10 also Ramachandran plots of alanine dipeptide are shown. The plots form the simulation should look like the right plot, but they look like the left plot. This can be explained by the force field used.

The CHARMM36 force field is an updated version of the CHARM22 force field and includes an addition called CMAP. The CMAP is an energy correction that improves the conformational properties of the peptide backbones.[1]

The Ramachandran plots in Figure 10 show the difference between the two for alanine dipeptide. The left and right plots are with and without CMAP respectively. As mentioned before, the plot without CMAP correspond to the plots in Figure 9. Thus, the alanine dipeptide simulations did not use the CMAP of the force fields.

Why this happened is not exactly clear, but it has probably to do with the naming of the atoms in the alanine dipeptide files. Because of the way the alanine dipeptide was made using PyMol, the names of the atoms had to be changed to correspond to the CHARMM force field. It is possible that the CMAP did not recognize these names and was therefore not used.

(20)

Figure 10: Ramachandran plots of alanine dipeptide. The right plot is without CMAP and the left plot is with the CMAP

6.2 3KI

Also from the PLUMED simulation on the 3K(I) peptide, Ramachandran plots were made. The plots are shown in Figure 11. Figure 11a en 11b, show the lowest and highest temperature replicas of simulations using the PLUMED replica exchange.

Figure 11c and 11d show the results of a REMD simulation done by my daily supervisor, Sander. Again the lowest and highest replicas are shown. The shape of the plots correspond to the plots with C-MAP in Figure 10, so here the CMAP is used in the simulation.

The high temperature plots look very much alike, so the same conformational space is explored. For the low temperature simulation, the REMD simulation samples more than the PLUMED simulation. An explanation for this is given in the following section.

(21)

-150 -100 -50 0 50 100 150 -150

-100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(a) Simulation 3KI PLUMED 275K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(b) Simulation 3KI PLUMED 427K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(c) Sander’s REM simulation at 275K

-150 -100 -50 0 50 100 150

Psi

<1 1-10 10-20 20-30 30-40 40-50 50-60 60-70 70-80

(d) Sander’s REM simulation at 427K

Figure 11: Ramachandran plots 3KI simulations

In Sander’s research[5], the melting point of 3K(I) with different water models was investigated. In Figure 12, the REMD results, PLUMED results and results from experiment are shown. Two things stand out in the Figure. One is that the helicity of the PLUMED simulation is higher then the helicity of the REMD simulation and even higher then the experimental data. Secondly, the decline of the helicity is much faster for the PLUMED simulation than for the other graphs.

The fast decline could be explained by the the way the PLUMED increases the temperature. The temperature of PLUMED is an ”effective” temperature and an ”effective”

temperatures of, e.g. 350K, does not necessarily correspond to the real temperature of 350K. So it is possible that the point at 300K should be at 350K. Only the temperature of the lowest replica corresponds to the real temperature.

That the higher replicas do not correspond to the temperatures does not make the simulations invalid. The higher replicas were created for more sampling, not to correspond exactly to the temperature. As long as configurations in the lowest energy replica is correct, the method is valid.

An explanation for why the PLUMED helicity is higher at 275K then the REMD helicity, is that the pressures of the simulations are different. The REMD simulation was

(22)

equilibrated with a pressure of 1 bar at 300K. To create a replica at 275K, a NVT run at 275K was done, using the result the NPT run at 300K. This would lead to a replica with a lower pressure at 275K. The PLUMED simulation was equilibrated at a pressure around 1 bar at 275K. A lower pressure leads to less helicity, since a lower pressure results in a less folded conformation. Thus the REMD simulation would have a lower helical fraction due to the lower pressure.

A simulation with less helical fraction can sample more configurations. This explains why in Figure 11c more space is sampled then in 11a.

Figure 12: Melting curve 3K(I). The black line is from experiments. the blue line corresponds to the REMD simulation[5] and the red line to the PLUMED simulation.

6.3 Hemagglutinin

During the HA simulations several configurations were seen. These configurations can be ordered in three groups. The different groups are depicted in Figure 13.

In 13a the starting configuration is shown, which was a straight, helical structure. In 13b, the peptide is in a folded configuration with two helical arms and in 13c a folded structure without alpha helices is shown. These configurations will be referred to as configurations A, B and C. The configuration shown in 13d, is from a crystal structure experiment[10].

(23)

(a) Configuration A, a straight helical structure

(b) Configuration B, a folded helical structure

(c) Configuration C, a folded non- helical structure

(d) Bend configuration from crystal structure [10]

Figure 13: a, b and c show different configurations seen during the HA simulation. d shows a configuration from a crystal structure.[?] The colors indicate the different types of the residues. Red indicates an acid residue, blue a basic, green a polar and white a nonpolar residue. In the Figure, histidine 106 and aspartate 112 are explicitly shown.

In Figure 14 the results are shown of the angle and d_{RM S} calculations, which were described in section 5.5. Also the course of the helical fraction is depicted in the Figure.

From these results and by analyzing the simulation in VMD, the course of the simulation can be determined. The course was as follows,

The simulation started in its straight, helical configuration. After about 150 ns the peptide began to fold. From this moment on, the peptide switched mainly between configuration A and B. Later, around 400ns, the peptide was always folded and started

(24)

to switch between configuration B and C. Until 550 ns, configuration B was most often seen, but after 550 ns configuration C appeared more often.

An important point about the simulation is, that it did not converged yet. In the last 100 ns the simulation is mostly into the non helical configuration, but the simulation would have to run longer to see, whether it stays in this configuration or goes back to a helical form.

(a) Angle calculations (b) Helical fraction over time

(c) dRM S per second (d) Distribution of the dRM S

Figure 14: HA analysis

Now will have a look at the d_{RM S} graphs. In Figure 14d, the distribution of the d_{RM S} is shown. It shows a clear peak at 0.3 nm. In a simulation done by Sander [4], starting in a folded configuration, the d_{RM S} was about 0.25nm. So, for an almost ideally folded simulation this would be the d_{RM S}.

The simulation of Sander was done with three whole HA monomers. Since the other monomers are excluded and only the hinge region of one monomer is simulated in the PLUMED simulation, some interactions are missing and a small deviation in the d_{RM S} is not strange. So the conclusion is, that the configuration B, comes close to the folded configuration from the crystallography experiment.

The results clearly show that the peptide goes into a bend configuration, even without the protonation of the histidine residue. This is in contrast with the results of Kalani et al.[13]. An explanation for the difference is the force field used.

(25)

As described earlier, the force field that Kalani and al. used, CHARMM27, has a tendency to be too helical. A more helical structure makes it more difficult for the peptide to fold, since the middle part must be non helical in order to fold. So, their neutral peptides did not go into a folded configuration, or even stayed folded, because they were too helical. In the case of the protonated simulations, the charged histidine disrupted the αhelices and the pepitde went into a folded configuration.

But why does the neutral peptide folds? When one looks at the configurations in Figure 13b and 13d, the nonpolar residues are on the inside and the acid, basic and polar residues are on the outside. So the hydrophilic residues are more in contact with water, than the hydrophobic residues. This leads to the conclusion that the bending is probably mediated by hydrophilic and hydrophobic interactions, not by the protonation of the histidine.

Another difference between my simulation and the Kalani simulation is the sampling of a non helical configuration. As said before, the peptide begins to become non helical after about 550 ns. The is not seen in the crystal structure or the Kalani simulations.

An explanation for the unfolding is that only a small part of the HA is simulated. This leads to a more flexible structure, which can lead to more unfolding.

The non helical configuration is not seen in the Kalani simulations, the reason for this, is probably that the force field they used leads to too helical structures, as is explained earlier. An other possible explanation, but less likely, is that no enhanced MD was used in the Kalani simulation. So maybe an energy barrier was overcome in my simulation in order to becomes helical.

The bending of the neutral peptide does not mean that the protonation of the histidine plays no role in the bending of the hinge region. A protonated histidine can disrupt the αhelical structure and make it easier for the peptide to fold. It would be interesting to see, how much faster a protonated hing peptide would fold in comparison with a neutral hing peptide.

7 Conclusion

To sample the large energy space for large molecules enhanced MD is still very important.

The first part of this study was about two enhanced MD methods, which made use of an implementation of Hamiltonian replica exchange.

The first method tried was the REST implementation suggested by Terakawa et al..

REST was tested on an alanine dipeptide and a 3K(I) simulation. Even though less replicas could be used, than with REMD, the simulations had a very large speed loss.

The loss in speed is largely due to the use of the computationally expensive free energy perturbation of GROMACS. The REST implementation is also very time consuming to implement without scripting.

The implementation of REST via PLUMED is promising. The speed of the simulations is significantly higher than that of the Terakawa implementation. For a 3K(I) simulation, the sampled configurations corresponds to the sampling of a REMD simulation. Also the helicity of the simulation comes very close to values from experiments.

Compared to REMD experiments. The method is very easy to implement due to avail- able scripts. The easy implementation and faster sampling makes the PLUMED method

(26)

favorable over the Terakawa implementation.

Second topic of the study was the bending of the hinge region of hemagglutinin. Even though the process is very important for the infection of viruses, the molecular details remain elusive. In a paper by Kalani et al. a mechanism is proposed and simulated where protonated histidines lead to the folding of the peptide, while neutral peptides remain unfolded.

In my study a neutral peptide was simulated using the PLUMED implementation of Hamiltonian replica exchange. In contrast to the Kalani paper, the peptide, starting in a straight configuration, did end in a bend configuration. The difference is probably due to the used forcefield. In order to fold, a part of the hinge must be non helical. The forcefield used by Kalani and al. was too helical, resulting into the straight, helical end configurations that were found.

So what makes the peptide fold? The folded structure from crystallography experiments and from my study show, that most non polar residues are on the inside of the structure and that the acid, basic and polar residues are on the outside of the structure in a folded configuration. This leads to the conclusion that the bending is probably mediated by hydrophilic and hydrophobic interactions, and not by the protonation of the histidine.

Unfortunately, the configuration did not yet converge to a steady configuration. For further research, the length of the simulation could be extended. Also a protonated peptide should be tested, to see what the influence of the protonation is on the bending.

Furthermore, could different mutations of HA be simulated. Simulations could also be done with the two other monomers present, so that the intra monomer interactions are also taken into account. And lastly, the simulation could be done with a larger part of the peptide, to get a more complete picture of interactions playing a role in the bending of the peptide.

8 Acknowledgments

I would like to thanks Erik van der Giessen and Patrick Onck for the supervision of my research and the helpful meetings. I would also like to thank Dogan Yilmaz for the installation of the PLUMED patch and for the solving of computer related problems. I especially would like to thank Sander Boonstra for his supervision of the project, the installation of PLUMED and for the help and guidance he offered daily.

9 Bibliography References

[1] Charmm tutorial. https://www.charmmtutorial.org/index.php/The_Energy_

Function", July 2016.

[2] M. Abraham. An introduction to replica exchange simulations: Mark abraham, session 1b. http://www.gromacs.org/Documentation/Tutorials/GROMACS_

(27)

USA_Workshop_and_Conference_2013/An_introduction_to_replica_exchange_

simulations%3A_Mark_Abraham,_Session_1B, July 2016.

[3] R. B. Best, X. Zhu, J. Shim, and al. Optimization of the additive charmm all- atom protein force field targeting improved sampling of the backbone phi, psi and side-chain chi(1) and chi(2) dihedral angles. Journal of Chemical Theory and Com- putation, 8:3257–3273, 2012.

[4] S. Boonstra. 3. 2015-07-10 fusion meeting - only chapter 2, November 2015.

[5] S. Boonstra, P. Onck, and E. van der Giessen. Charmm tip3p water model suppresses peptide folding by solvating the unfolded state. The Journal Of Physical Chemistry, 2016.

[6] B. Brooks, C. Brooks, A. Mackerell, and al. Charmm: The biomolecular simulation program. Journal of Computational Chemistry, 30:1545–1614, 2009.

[7] B. Brooks, R. Bruccoleri, B. Olafson, D. States, S. Swaminathan, and M. Karplus.

Charmm - a program for macromolecular energy, minimization and dynamics calculations. Journal of Computational Chemistry, 4:187–217, 1983.

[8] G. Bussi. Hamiltonian replica-exchange in gromacs: a flexible implementation.

Molecular Physics, 112:379–384, 2014.

[9] G. Bussi. Hrex in plumed. https://github.com/GiovanniBussi/plumed2/blob/

v2.1-hrex/user-doc/tutorials/hrex.txt, July 2016.

[10] J. Chen, J. Skehel, and D. Wiley. N- and c-terminal residues combine in the fusion- ph influenza hemagglutinin ha(2) subunit to form an n cap that terminates the triple-stranded coiled coil. PNAS, 96:8967–8972, 1999.

[11] U. Ess, amm, L. Perera, M. Berkowitz, T. Darden, H. Lee, and L. Pedersen. A smooth particle mesh ewald method. 103:8577–8593, 1995.

[12] S. C. Harrison. Viral membrane fusion. Journal of Virology, 479:498–507, 2015.

[13] M. R. Kalani, A. Moradi, M. Moradi, and E. Tajkorshud. Characterzing a histine switch control ph-dependent conformational changes of the influenza virus hemagglutinin. Biophysical Journal, 105:993–1003, 2013.

[14] C. M. Mair, T. Meyer, K. Schneider, Q. Huang, M. Veit, and A. Herrmann. A histidine residue of the influenza virus hemagglutinin controls the ph dependence of the conformational change mediating membrane fusion. Journal of Virology, 88:13189–

13200, 2014.

[15] S. Marqusee, V. H. Robbins, and R. L. Baldwin. Unusually stable helix formation in short alanine-based peptides. Proceedings of the National Academy of Sciences of the U.S.A., 86:52865290, 1989.

[16] R. Nifos and F. Oteri. Cmap addiditon partial tempering script. https://groups.

(28)

[17] P. Palese and M. Shaw. Orthomyxoviridae: The viruses and their replica- tion. http://fly.reactome.org/cgi-bin/eventbrowser?DB=gk_current&FOCUS_

SPECIES=Homo%20sapiens&ID=168275&", June 2016.

[18] Schr¨odinger, LLC. The PyMOL molecular graphics system, version 1.8. November 2015.

[19] T. Splettstoesser. Landscape. "http://www.softsimu.net/Talks/Pics/

landscape.png, June 2016.

[20] T. Terakawa, T. Kamedi, and S. Takada. On easy implementation of a variant of replica exchange with solute tempering in gromacs. Journal of Computational Chemistry, 32:1228–1234, 2011.

[21] G. Tribello, M. Bonomi, D. Branduardi, C. Camilloni, and G. Bussi. Plumed2: New feathers for an old bird. Computational Physics Communications, 185:604–2013, 2014.

[22] G. A. Tribello, M. Bonomi, D. Brandyardu, C. Camilloni, and G. Bussi. A quick introduction to plumed 2. https://www.youtube.com/watch?v=PxJP16qNCYs, June 2016.

[23] D. van der Spoel, E. Lindahl, B. Hess, and the GROMACS development team.

GROMACS User Manual version 4.6.3. www.gromacs.org, 2013.

[24] L. Wand, R. A. Friesner, and B. J. Berne. Replica exchange with solute scaling: A more efficient version of replica exchange with solute tempering (rest2). Journal of Physical Chemistry, 115:9431–9438, 2011.

(29)

10 Appendix

How to apply REST in PLUMED

The first part of performing REST in PLUMED is to create the different replicas that will be used for the replica exchange. This section is largely based on a tutorial written by Bussi.[9]

To create different replicas, the so called ”partial tempering” script can used.[21]

This script is a part of the PLUMED program and modifies the necessary force field parameters. The script makes use of a processed topology file, which is created from the topology of the lowest energy replica. A processed topology is a topology file, in which all the files that are included in the normal topology file are written out explicitly. So for example, all angles and dihedrals parameters are written out. Now the steps necessary to create the replicas are discussed.

- The processed topology file is created from the original topology file by using the grompp program in GROMACS with the -pp flag.

- In the [atoms] section of the processed topology file, an underscore, ” ”, must be added to the atom types of all atoms in the solute.

- With the processed topology, including the underscores, the partial tempering script can be used. The script creates a processed topology file, in which the necessary parameters are scaled with a factor $i. The following syntax is used to create replica

$i,

plumed partial_tempering $i <processed_.top> topol_$i.top.

- The script is not entirely complete. It does not scale the CMAP and the pair types, which also need to be scaled for a correct implementation of REST. In the following syntax the CMAP is also scaled.[16]

plumed partial_tempering $i <processed_.top|

awk -v f=$i ’ BEGIN {cmaptypes=0; pairtypes=0}

{

if (/\[ /) {cmaptypes=0; pairtypes=0};

if (/\[ pairtypes \]/) {pairtypes=1};

if (/\[ cmaptypes \]/) {cmaptypes=1};

if (pairtypes==1) { print $0;

if (substr($1,0,1)!=";" && substr($1,0,1)!="[" && NF>0)

printf "%s_\t%s_\t%d\t%.12lf %.12lf ; scaled \n",$1,$2,$3,$4,$5*f;}

else if (cmaptypes==1 && substr($1,0,1)!=";") { sub(/\\/, "", $0);

if (NF==8) printf "%s %s %s %s %s %s %s %s\\\n",$1,$2,$3,$4,$5,$6,$7,$8;

else if (NF==10) printf "%.9lf %.9lf %.9lf %.9lf %.9lf %.9lf

(30)

$4*f,$5*f,$6*f,$7*f,$8*f,$9*f,$10*f;

else if (NF==6) printf "%.9lf %.9lf %.9lf %.9lf %.9lf %.9lf\n",

$1*f,$2*f,$3*f,$4*f,$5*f,$6*f;

else print $0;}

else {print $0;}

}

’ > topol_$i.top.

- The script is not very robust and errors can be made, so it is advisable to check if the scaling was done properly. To do this one can do two tests:

* Run the ”partial tempering” script with a scaling of 1. Then do a rerun of a trajectory with this scaled topology and an unscaled topology. The results of the simulations should be the same.

* The same can be done using a scaling of 0.5 instead of 1. Now the relevant energies, like the Lennard Jones and the dihedral energies, should be half for the scaled simulation.

Now that the replicas are created, the HREX can be performed. The process is not much different than for normal REMD in GROMACS. A difference is that a plumed.dat file must be included. In this file one can setup different options for PLUMED. With HREX this file is not used and an empty plumed.dat file in the map with the replica is sufficient. The following syntax can be used to perform HREX in PLUMED, with the replicas in different maps,

mpirun -np 30 --output-filename mdrun.job mdrun_mpi -v -plumed plumed.dat -s topol.tpr -multidir Map_1/ Map_2/ Map_3/

-replex 500 -hrex -cpnum -cpt 180 &>mdrun.job &.

The bending of hemagglutin