 Hi everyone, my name is Nora Aho and I come from the University of Jyväskylä located in Central Finland. And today I want to talk to you about something called constant pH molecular dynamics, that I am currently as part of my PhD implementing into the Gromax software. So I am working in the research group of Gerrit Grönhoff together with Pavelupus Live and we are closely collaborating with the Gromax developing team in Stockholm. So in KTH Stockholm we are working together with Berkes, Paul Bauer and also Anton Janssen to get this implementation done. My research funding comes from the Academy of Finland and also I am using the computational resources by the CSC computing center in Finland. First of all, why do we need constant pH molecular dynamics? So I want to introduce you briefly a background of why is this important. So maybe you know that pH is one of the parameters, one of the main parameters that affects the dynamics in biological systems. So when pH of the environment of the solution changes, this leads to changes in the protonation states of certain titratable groups and also that then changes the charges of these groups. Change in protonation, for example in proteins can be directly coupled to structural changes. So when the charges change, all the interactions change and even the whole structure of a protein can change. Still, pH is not something that one can easily take into account in molecular dynamic simulations. Unlike pressure and temperature, pH is not a parameter that one normally fixes to be constant. And this is where constant pH MD comes into hand. So it is one of the methods that then would enable us to take pH into account when running molecular dynamic simulations. In my talk, I want to briefly introduce you to the method that we are using that allows one to run constant pH molecular dynamics. And also I want to give you a kind of user point of view on what would you then need to take into account if you after this talk or in the future feel like pH is something that could be important for your specific system. Okay, so the way we have been implementing constant pH molecular dynamics into gromax is with a so-called lambda dynamics approach. This approach was first introduced by Charles Brooks already 25 years ago. And the idea is that if one has a titratable group here, the picture is from glutamic acid, then one can assign a coordinate called lambda to describe the protonation state of this group. So here, for example, we set lambda zero to describe the system that is protonated, the hydrogen is there. And as lambda then goes to one, the charges would be deprotonated and the group would be deprotonated. So the idea of lambda dynamics is to connect two different thermodynamic states via a lambda coordinate. And here, unlike, for example, in free energy calculations, this lambda is not fixed, so it will not have fixed values in this interval, but instead it is a continuous variable that we also assign a mass. So it can freely move between these two states depending on the interactions with the environment. One thing to note here is that this hydrogen at lambda is one when the group is deprotonated in our implementation, it's really not getting off. So the hydrogen atom will physically still be here, but only the charges of the hydrogen will be zero here. And the other charges will be the deprotonated charges. So we are not really breaking this bond, but just the charge interactions are gone. Okay, so then what these lambda dynamics and the constant pH adds to the normal standard MD is that each simulation step, we also need to update the values of these lambda coordinates. And that then changes the charges of tight-radable groups. The way we start thinking of how we add the new lambda degrees of freedom to the total Hamiltonian of the system is that we add a few terms to the total Hamiltonian. So instead of only taking into account the Cartesian coordinates, now we also write the Hamiltonian as a function of lambda, because when lambda changes, the charges of the system change and this changes the interactions. So for constant pH Hamiltonian looks something like this. There is the kinetic energies of both the new lambda particles and the Cartesian coordinates. Then we add a few terms that only depend on lambda. These are just specific for constant pH MD. There is mainly, the most important is the effect of pH. It is a linear term that depends on the reference pKa of your group. So when pH is set to a value, it can then tilt the energy landscape and start favoring one of these protonation states. We also need to add a free energy of the protonation term just because we need to describe the force field term of breaking this bond. Then we add a third term, which is called the biasing potential, which just allows us to sample from the interval lambda 0 to 1. So we add high walls to prevent sampling other values of lambda. These are all quite easy to handle, but the most, let's say, challenging part is the potential energy of the system, because this now also depends on the lambda values. So let's look a bit more into details of that. The potential energy of your system consists, as in normal MD, of bonded interactions and non-bonded interactions. So there is then van der Waals and Coulombic interactions. In our implementation of constant pH, we only take into account that when lambda changes, the charges change. So we do not consider bonded interactions or van der Waals interactions and we are not considering that the atom types in your simulation system are changing. This is also an approach that was used before and is used by other groups. So we only consider the Coulombic interactions for the charges that are also changing. And the Coulombic interactions for the atoms that can change the charge for the atoms that are part of a titradable group is described by interpolating between the two states. So let's look a bit more details into that. So we describe the Coulombic interactions by interpolating between protonated and deprotonated states for the groups that can be titradable. And the Coulombic interactions, the Coulombic potential looks something like this. So if one has one lambda group, so you consider that in your normal simulation system you add one group that can be titradable. Then you will have some different terms in the Coulombic potential. So there is a term that looks like the standard Coulombic interactions where you have all the atoms that have fixed charges. Then you add two terms out of which the first one considers Coulombic interactions between the fixed charges charged atoms and then atoms that have that are titradable. So all the atoms that are titradable in the Coulombic interactions will then have a charge that is described by a linear interpolation between Qa and Qb. The charge is for protonated and deprotonated group. And then the term which describes Coulombic interactions between the atoms that are both part of a titradable group will look something like this. So we'll have two of these interpolated charges. And what we then need to calculate for each step of the simulation, we of course want to calculate the force acting on the lambda particles because that is how we can also update the lambdas. So in addition to calculating the force acting on the atoms, we also calculate the forces acting on these lambda particles. And this can be calculated by using the electrostatic potential and the charges of each atom. So when we are running constant pH with our implementation, we do not need to add any other calculations on top of the standard force calculations. We only need to extract with the forces the electrostatic potential when we are calculating these Coulombic interactions, for example with PME. And this is something that is really important because in the previous implementations, so there exists an implementation of constant pH MD for chromax from maybe 10 to 5, 5 to 10 years ago. But there the problem was that these interpolation terms were done in a slightly different way, which led to multiple force calculations if you added more and more lambda groups. But now when we are using the electrostatic potential, the only thing that changes is that we need to extract this. So we still only need one force calculation as in standard molecular dynamics. And that I am showing here in this graph, because one of the main things of our new implementation of constant pH is that it can in principle be as fast as standard molecular dynamics. So previously, in the old implementations for gromax, the issue was that the more titratable groups one adds, the less the smaller the so the simulation speed was decreasing. But now as you can see, even if you add here 16 or 100 titratable sites, the speed of your simulation does not go down. And this of course allows us to sample bigger and bigger systems. And we are really happy that there is this new way of interpolating and writing these equations because that allows us to get a really nice speed for also the constant pH simulations. One additional thing, because in the previous examples I was only considering a group with two different states, but of course there are titratable groups that have multiple protonation states. So in our implementation, together with the speed up, together with the possibility to run as fast as standard MD, we also introduce a new scheme to treat the multi-state titratable group such as he's been shown here. And you don't need to know the details right now, but if you are then interested, you can later look into details. So the idea is that one assigns one lambda coordinate per protonation state and then these lambda coordinates are coupled together via a constraint. The constraint is similar to the shake constraint used for treating the other constraints in your simulation system. And then one is able to reproduce, for example, for histidine the titration curves. An important to notice is that having these multi-state groups also does not decrease your simulation speed. So in our current implementation it doesn't matter how many groups one has and how many protonation states these groups have, you can still in principle run as fast as standard MD. Okay, so in practice I want to briefly tell you a user point of view story of how one would then set up constant pH. So prior the simulation the user needs to select the interesting groups to treat as titratable. So for example you might have a few groups in your protein that are, that you know from experiments that are interesting and you would select those. You need a reference pKa for each of these titratable groups and charges from the force field. That means that for example for amino acids of course you know the reference pKa values and we are also providing this information, but if you have a more specific group maybe you need to find out yourself what is the reference pKa. Another thing that one needs, you need to know this term of free energy of deprotonation. This is something that you can compute similar to the normal free energy calculations and this needs to also be computed for a reference group in water. So for amino acids we can provide these parameters for the free energy of deprotonation for standard groups for amino acids in water, but if your system consists of something more special then you need to think a bit on how to run these calibration runs and get these parameters. Then the next thing is to select pH, some more or less default parameters and then run at your selected pH. So a simple example would then be to take for example in the cardio toxin protein, take four different groups that you think are interesting, set them as titratable, you assign one lambda for each group, okay for this histinine you assign a bit more lambdas and then you run constant pH simulations starting from the same structure at different pH values. Then similarly to titration experiments, yes titration experiments you can then from your simulations extract the fraction how much your group would be during the simulation protonated or deprotonated. So you can calculate the fraction for the group favoring to be in one of these protonation states and from this information you can fit the titration curves to your data. In addition to getting pKa values from your simulations, of course this is a nice way to validate the constant pH method and it's also a nice tool to get this pKa information, but in addition one can then also study the conformational changes in the protein. So there might be a case when where a titratable group changes the charge, changes the protonation state and then the whole structure of the protein is changing and this is the kind of information we are also hoping to use the constant pH4 in the future. Last I want to tell you the real status of our implementation because maybe it already sounds good that okay let's start using it, but in reality we are not quite yet there. So we have, I have implemented an improved constant pH scheme to the main code of Gromax with the lambda dynamics-based efficient method, but currently we are in this testing phase where we still have a few challenges that we are now trying to fix. So for example we still have some challenges related to the multi-state groups and also we are coming up with their routine to keep the total charge neutral because this is something that is important, the more titratable groups you add that your, that the total charge of your box is not fluctuating and what we are still yet to do is to really publish this and put the code into the official Gromax release. It's something we are still hoping to do this year and then in the end release the code in the official Gromax together with the user guide, tutorial and also example systems. Okay, I want to tell you last two things maybe to remember from my talk is that constant pH MD is a method which allows one to include the effects of pH into simulations and in addition to configurational space one is also then sampling the protonation states and the second thing is that constant pH MD can be as fast as standard MD and that is important for running bigger systems. I want to thank you all for listening and you can go ahead with your questions then live with me in the Q&A session or you can also always send me an email and the address is shown there. Hope you enjoyed.