 Thank you for inviting us. Thanks, John, for having us over, who's not in the room yet. But I'm very pleased, and Jeroen and I were both very pleased to be here and talk to you about some of the things that we have thought about in terms of force field development for the Charm Force field. And a lot of the ideas that I am going to talk about, maybe old news for you, because a lot of the stuff that we use and that we thought about are inspired by work that has been done in the context of the open force field initiative. So before I start, I want to give a brief outline. First of all, I want to talk about the Charm Force field, the Charm Lipid Force, what it is, why we're interested in it, and what we want to do with it. That's going to be the part of our dramatization. And then I'm going to present to you two things that we've already done. One is the improvement or improving the description of permeabilities in the force field. The other one is incorporating long range and Jones interactions. And I think or I hope that I'm going to be able to give you a taste about lipid force fields and why they are hard to do, like why it's hard to do the good lipid force field. So of course lipids are the basic building blocks of biological membranes. And in order to simulate lipids, we have to describe a lot of interactions between the lipids, between lipids and water, lipids proteins, lipids and small molecules. And in order to do this, we of course want to have a force field where we get all these force constants and all these constants in the force field, the force field parameters, so that we have a reasonable description of the atomistic interactions. And the Charm Lipid Force field has been the most used lipid force field over the last 10 years. And since its inception in 2010 by Jeff Claude, it has according to Google scholar over 2000 citations. And when we think a little bit about the reasons for its success, then I think one thing is that the Charm Force field is a pretty general solution. So a pretty complete solution. There is a protein and a DNA lipid and they all are pretty good. And then there's the Charm Generalized Force field for small molecules. I think the second reason is Charm GUI. This input generator can generate membrane proteins embedded in membranes and all this stuff online, and not only for Charm but also for other programs. And that's why I think the Charm Force fields are not only popular in the Charm community, but also in people in Gromax and Amber, they use the Charm Force field just because it's convenient. I think a third reason is probably its stability. That's mostly thanks to the consistency and thoughtfulness of Alex McCrell, not letting everything into the force field, but really having kind of a steady improvement. And the fourth reason for the Charm Lipid Force field being popular is that it's pretty good. Like it does a very good job in getting many important properties right. And those are mostly equilibrium properties of lipid bilayers. So when you, or biological or bilayers and biology are self-assembling, so which means that they have an overall zero surface tension. And when you look at the the tangential pressure profile in a membrane, there's two lipids that with the tails pointing to watches at each other, then you have that huge negative interfacial tension in the water head group interface, and then a huge tangential pressure in the chain, like in the midplane. And these two contributions to the tangential pressure have to zero. They have to vanish in order for the membrane to be self-assembling. And this is, this is, of course, depends on all the little interactions inside the membrane. It's really hard to get. But you have to get this right in order to get, for example, the surface area of the membrane brain right to get the right compressibility moduli to get good bending constants and spontaneous curvature. And the channel force field does all that, the channel bit force field does all that pretty well. When we look here, the simulated versus experimental area lipids for a lot of important lipids, simulate that as very close to experiment. And same for the bending constants. So C36 is great. A lot of people like it. Why should we bother optimizing it? I want to talk a little bit about what it doesn't do so good. The first thing is that interfacial systems in general are better described by long range than a Jones ejection. So dispersion interactions have an important long-range effect and interfacial setups. And like a few years ago, or for the last couple of years, most important simulation programs have an efficient LJPME implementation where you can use long-range and Jones forces and it doesn't really hurt your simulation speed. And the C36 force field for proteins has been validated for LJPME. So the protein force field works, but the bilayer force field, because the bilayer is so sensitive to all these interactions, when you apply it to a bilayer, you see that the surface area goes down. It's consistently down because for different cutoff distances. So what we're hoping to get from LJPME is that this line is constant, which it is, like the surface area does not depend on the cutoff radius. That's what this long-range LJPME should give us. And when we have a usual cutoff scheme, then the surface area decreases towards that limit. But the force field has been optimized for a cutoff distance of 12 angstrom. And you see that this is pretty good, but when we take all these long-range dispersion forces into account, then the force field doesn't perform as good for the area for lipid and for many other important properties. So the cutoff radius is an essential parameter of the force field. And we want to get rid of that. We're going to get all this long-range magic into the force field. One more thing that is not so great is that the bilayer areas, bilayer surface areas are good, but the monolayer surface tensions are really, really terrible in the KME force field, unless you use a long-range LJPME. Then suddenly, the monolayer areas are good, but the bilayer areas are not good more. So there is this thing, like an inconsistency in the force field. And there's also inconsistencies regarding Leonard Jones' treatment. For example, the protein force field was parameterized with a switching function between 10 and 12, the bilayer between 8 and 12, and some lipids in the bilayer in the lipid force field also from 10 to 12. So there are some inconsistencies that would all go away if we get rid of all these cutoff schemes. So if we just say, okay, we don't want to have to do with Leonard Jones' switching and with cutoff radii, we want to replace that by something that gives us long-range dispersion. Then the second thing that's not good about, or maybe not so good about the Charm force field was tip 3P water, which we all know is a crappy water model, but it's been in there from the start and we're not going to get it out. But it's really like almost every compound in the Charm force field is optimized for tip 3P. And replacing it by a better water model is going to be a major pain. And then the third thing that we found at C-36 doesn't do so well as permeability. I don't know how many of you have worked with in the context of the sample competition. But I mean, permeability is basically a partitioning plus some kinetics. And like with these MM force fields, we don't even get the partitioning right. So I mean, the success or the performance of MM based approaches in these competitions is usually pretty cool. So permeability is another thing that we want to optimize that the current Charm force field doesn't really get and no MM force field really. So next I want to talk about how we approach this these problems. And this is very much intertwined with my personal history or with my personal journey in science. I did undergraduate research with in a group and journey that was mostly doing phosphate parameterization using automatic approaches back in 2011. And yeah, this is me back then. And then I went on to do a graduate research in computational fluid dynamics. So something completely different. But I want to talk a little bit how this guy viewed the world and what he thought about automatic parameterization. So the consensus in our group back then was that, okay, a lot of people do force field development by hand. And this in principle is something that is very much suited to automatic procedures. So we were thinking about how can you use numerical optimization algorithms to improve the force field development process. And this would like the hope is that when we do this automatically, then we can make the development much faster. We can also cancel out human error. And we can make the force field more maintainable. So when you have the force or when you have a human design the force field, then it's really hard after 10 years to go back and improve it. Like for example, tip three p it's in charm. And it's going to be in charm indefinitely. And when you think about this, those are a lot of the problems that we deal with in code development. So in code development you also have this potential of human error where you make like you introduce one back and everything or nothing works anymore. Or you where you have like a huge piece of software that needs to be maintained over many years and the software development industry has found ways around this and ways to deal with these problems. And the hope is that when we automate the force field development process that we can like what can we learn from the software development industry. So what I mean basically what we used to do is we would start off with a set of force field parameters, run a simulation as a black box, get the observables, compare them with some reference data and do an optimization step to get a new promising set of force field parameters. So that's like 10 years ago or nine years ago. And then we the group that I was in used mostly gradient based optimization algorithms that we're operating on this black box simulation. And of course when you do this you have to like for each gradient calculation you have to run a simulation and given that the simulation is a black box for us we have to run another simulation to get this finite difference step and to calculate a partial derivative. And just in order to get a gradient at a given per average set you have to put in so much work and what you're going to get out of it is not like it's going to be very noisy because when you think about like two or the statistical uncertainties on these on these simulation results when you calculate a gradient or a finite difference between those who knows where this gradient is going to take you. And what I did in my master's thesis was to replace the gradient based algorithms by a meta model based algorithm. We basically say we don't want to get this local information instead we predict the outcome of a simulation using everything that we know from previous evaluations or from previous simulations. Now this is still like viewing the simulation as a black box but we instead of having this local gradient we have like a global optimization to our simulated data and then we use this meta model this predictor to detect new promising parameters and after simulating these new parameters we have something that we have new information to enhance the meta model so like these two things the meta modeling and the optimization they complement each other. And then you need some kind of of sampling technique that says okay I have a given model of my simulation or a meta model of my simulation how am I going to use that to detect new promising parameters and this worked well for when you have very few parameters like maybe eight to ten parameters in your force field and it worked well for a few things that we did with small molecules but you're not going to get meaningful predictions for when you have a high dimensional parameter space. Okay so this is a like a small toy example we have some some loss functional then we just randomly sample a few initial parameters from very short simulations for these and then we can get like an idea about how our loss function looks like how our observables looks like depending on the parameters and we can use that to sample these more promising regions of the primary space and remember this like the philosophy here was simulation is a black box remember everything you can build this model and then you can do this faster than with gradient based methods and then when I started with Bernie and Rich in 2017 like two years ago I took some of these ideas with me but I was also learned learned a lot of new things and when I first read these like of course a lot of a lot had happened in the years in between and when I first read the force balance papers and the m bar free energy things that just blew my mind because like coming from that background where you have to evaluate gradients and these finite differences and then reading about hey you don't have to run a new simulations for this and the simulation suddenly is not a black box anymore that was that was really huge for me to understand so I mean the basic idea most of you will know this is that you have a simulation trajectory that you generated in some ensemble a which corresponds to our original parameter vector and then you perturb your parameter vector which gives you a new ensemble b and so what we did was run new simulations for all of these particular parameters but you can when you write this as a thermodynamic dynamic average and multiply this by like interpret this in the in the original ensemble then you get this reweighted expression where you can use the original trajectory to calculate these partial derivatives without any more simulations and this is dirty and accurate and as I said probably old news to most of you but for me this was huge when I when I read this two or three years ago and as a closing or as a summary of this automated optimization thought when we think about numerical optimization it's basically metamodelling so we can decide whether our metamodel is global or local so basically all when you think about Newton optimization or steepest descent or conjugate gradient you learn something about the local shape of your function learn the gradient and you use that to get your next iteration and this was what this this grow method this find a difference method and then we what I did was use a global metamodel and set and use that to do a numerical optimization then you can of course do this physics informed we use reweighting as a local metamodel which force balance does or you can get the same from multi-state benefit and then the question is how can we get a more global metamodel because these reweighting approaches only work for small regions of the primary space where you have good confirmation overlap how can you make more global predictions and michael's shirts group had this paper about this about the pcfr approach and one question is how like how can you maybe use information from your local gradients to construct a more global metamodel without enforcing confirmation overlap so kind of combining these two ideas okay so with that I want to talk about two things that we have done over the last two or three years and the one is improvement of permeability description as I said in the beginning primabilities are basically partition coefficients plus kinetics so we can express them through the free energy profile like a barrier or well through the membrane like this this is the membrane and a diffusion constant so how fast does a particle move at a given point in the membrane um and with with m m and protest we usually don't even get this free energy right and we can't even say something about the diffusion profiles that we get because don't think there is experimental data to support this um or to to get the same detail or same precision of information but yeah one thing that's important here is that the dominant contribution to permeability is this water oil partition coefficient and so we we thought about we thought that okay in order to get better permeabilities the most important thing to get right is this water oil partition coefficient and um this is something that is going to be really hard to get in a general framework so we thought about how can we use or how can we develop special um special solutions that work for given molecules or for given permeants um so one one example of this is water partitioning into membranes so we look for example at tip 3p going from the water phase into oil into hexadecane and we saw okay the free energy of transfer is overestimated by 1k kelvin mole okay tip 3p is a crappy water model um what happens if we use something better like tip 3p force buttons or tip 4p force buttons or opc and it turns out the better the water model gets the worse this transfer becomes um and that's because this like when when you when you optimize a water model you're explicitly optimizing it for the bulk phase and you're you're you're very intentionally putting this into the optimization that okay we're going to describe bulk we're not like you have usually have this gas phase dipole correction you're basically saying i don't want to get this right i'm not interested in water translocate or water translocating into um a low dielectric phase and so when you when you look at water permeability it's going to be too high because you don't even like you don't even get into the oil phase as much as you would like an experiment and um yeah so so we thought about what can we do about this if this is bad for water it's it's even worse for more complex molecules um and then oh yeah and then sure enough so the first thing is what happens if we include polarizability if we if we allow the water dipole moment to relax in the hydrophobic phase and polarizable model water models don't have like this this fitting or that they're supposed to work in different kinds of environments that that's what the polarizability is all about right that a molecule can respond to different electrostatic environments when we looked at the partition coefficients for um for the swim form swim six group water models then this is almost perfect so this is um within um experimental error um and the question is how can we get this or can we get a similar result in an additive framework or can we do about the additive force field and then we develop this very crude procedure where we say okay say we have solubility data from experiment or from cosmo methods or from anything that can do it better is better than mm and how can we use this solubility data um to tweak um our interactions and basically tweak um the or individually tweak the lana john's interactions between water uh and the solute and between alkane and the solute the thing is that that we can do this very very cheaply so usually when we when we calculate solubilities we have these chemical growth simulations where we have different lambda states and then when you when you think um like the three final lambda states give you a lot of different configurations for partially um partially partially annihilated versions of the water and alkane we can just use that as a basis to do reweighting um and to i'll be basically use these the the m bar predictions from these three lambda states to generate lana john's parameters between water and alkane like pair specific lana john's parameters that will um give us the correct free energy um and this and when you use these kind of m bar metamol as you can do this in like in a day on a gpu easily so i think this this took us five five hours on a single gpu this whole optimization procedure plus one final simulation um okay and then we like we did this for water and we said okay if it works for water because i mean water alkane or water lipid interactions um like when we tweak with when we when we play with the water alkane interactions then a lot of things can go wrong potentially because as i said the balaya is so um so sensitive to these things and if we can get water alkane right then we probably can get everything else right as well um in terms of the of the permeability and for these optimized parameters we then measure permeabilities um through um some exemplary balayas and you can see those are the this is the original c 36 phosphate with chips 3p and then when we use these optimized parameters then um it's a little bit or it's it's significantly closer to experiment but still not right this is this kind of a um a solution like an intermediate solution ultimately you want to do this in a physically more rigorous way um but this uh works and gives us the correct permeability for um for water and you can and i mean this is directly transferable to other molecules as well um i'm not doing any time yeah plenty of time okay all right okay so now i'm going to talk about this um adapting charn to long range dispersion how can we tweak the parameters of the charm force field to get all the things that i talked about in the beginning right like the this inconsistency between balay and monolayer and um how can we or the the incorporation of long range dispersion into the charm force field um and this and you're loondered a lot of this work um oh internet connections in i think it's fine yeah it's fine okay all right you're loondered a lot of this work um and we've been working together on this for a year or a little bit over a year um so the objective is we want to adapt the c 36 parameters to work with lj pme and we want to get get get the same quality as charm 36 for the properties that charm 36 does well which is aerial billet bid and the um the scd's the auto parameters have compressibility like all these structural properties for bilay is we want to retain but at the same time we want to achieve consistency between monolay is and bilay is um and we also don't want like in a in a lipid the tail parameters are basically alkan parameters um and then the head group are specific lipid parameters and we didn't want to touch the alkan because we've validated that like we um add a similar dnr um published a paper about this last year um where she calculates all kinds of properties and she consistently shows when you use lj pme you use long range jenna jones and all these parameters become better um yeah i'm not gonna go into this table yeah maybe maybe i should maybe i should this just shows how like this is from the original c 36 paper and this shows surface tangents for c 36 monolayers with the default um settings for the charm force field you can see that they are they are much lower than experiment and when you use like a long range dispersion and they become not perfect but much better um but this is the kind of inconsistency so so so these are standard charm settings these have non-standard settings and we want to have a parameter set that works for the standard settings or yeah for like for want to have a set of parameter that works with a consistent setting for both monolayers and bilayers okay so our optimization targets are the the area of both monolayers and bilayers and then the compressibility um because the compressibility is is something that is directly influenced by um by this long range dispersion and it's also something that is that can go wrong very easily as is seen for example in the charm um in the polarizable lipid force field that has pretty good areas but the compressibility is just not very good yet um our next observable is the auto parameter which describes the torsions of all the lipid atoms in a bilayer and then um we there is also experimental data for the hydration of the head crew but we don't like we would we're taking this as part of the optimization because it's cheap we don't have to run um a full bilayer or monolayer for that um but we're also not not waiting it so much um because like we know that tip through p is not going to do a super good job in getting all these hydration things right all the details of the hydration right okay and then we can use um reweighting to calculate the sensitivities um of in this case the area prolipid towards um certain parameters in the simulation so for for the charges um we see that or so the first thing that we see is that um it only works in a very small range for a complex system like that and that's because when we move out of this range we don't have enough confirmation overlap and we'll just get basically a constant and very large error bars so these these error bars are always over three replicas of the simulation um we also see that for charges and um yeah like for all these charges we can change the charge quite a bit like 0.01 still get a reasonable prediction and it stays kind of consistent whereas for um a lemma jones sigma it's it's much harder because conformational overlap or conformational overlap is more of an issue and um for the lemma jones well depth um you it's it's it's almost linear over a certain range if the sensitivities are almost linear over a certain range and the um the prediction uncertainties are very consistent predictions are very consistent between the replicas because because you're not messing with the distance you're just messing with with the energy of the interactions and then but then the question is how well do these m bar or reweighting predictions um in comparison to real simulation so what if we simulated with um or if the change the parameters are simulated and we see that in in a small range around the original parameters the simulation and the reweighting agree very well um for for charges and they at least give the same trend for sigmas so one thing that's specific for liquid force fields is that when you change change the lemma jones radius in this case i think this is the cardinal oxygen either this is the no this is the phosphate oxygen um and you when you change it only by four percent or so then you get a phase transition like this is one parameter in the simulation and like the area depends like first of all what's what's surprising is you increase sigma and and the area like you increase the radius and the area of the thin strings because there is no more not enough hydration between the head groups because now now the the liner jones interacts are taking over um between the lipids not so much with the water that's hydrating that root and then when you change it only like by five percent like one parameter is in the simulation the area goes down by ten percent and you get a phase transition in terms of like the from the liquid crystalline phase which is very pretty mobile to a gel phase that is that is more rigid and more dense and and something like that you can obviously not capture by a reweighting approach but also hear you at least get this like the direction of like you you get an idea of in which direction the parameter changes so do you see that the effective sample size plummets as you make that change so that you can detect that the weights are all going to zero except for one weight one snapshot so we've been trying to find automated ways of detecting when we're too far to be able to extrapolate and we often find that first of all the errors become big and then they become small again because you end up using just a few samples from your ensemble to read to capture the extrapolate property so if you uh do you find that the weights indicate that you have a very small effective sense so i don't think we use um far for these calculates we just use the one side reweighting because it's basically just one state right because you're computing an effective way to each of the mice over uh the you know normalizing constant and and we haven't looked at that but um the variance of those weights would go way up and then one of them would be uh very large the rest would be small what would be you can take a look of some of the masses that as you get too far away then um this is just a general feature of any sort of reweighting scheme the effective sample size plummets because just a couple of those that have the energies that are lowest at the new parameter set end up dominating by such a large energy gap so uh monitoring that it turns out to be a decent way but it's still heuristic you ever seen when when this divergence should happen one way to think about this is for all these reweighting methods you have uh you know each sample is associated with a weight so instead of being one it's the weight and you can convert some or this should already make us think and think like okay getting like these small details have such a large effect um that's that's kind of demonstrates why getting a good lipid force through this hard but um then you have also this this very non-linear effects for example um for this for this oxygen uh like for this oxygen charge here when you change the parameter or when you change the charge then uh the parameter doesn't or the the area of lipid doesn't follow linearly but it has some some more complex complex behavior behavior um and still the reweighting gives us a feasible decent a decent direction for the for an optimizer but the but the dependence between parameters and observables can be very very non-trivial and we see that here so when we change the charges by zero zero five electric charge elementary charge units then um um like a lot of different things happen like the area can can can change up to one angstrom squared um but um it's very hard to predict as a human what would happen and that's why i think at this level of fine tuning we really need these automated approaches a question can i ask a question quickly this is what i talked about earlier actually most line of jones radii in the head group when you decrease them or when you increase the the the radii then the area collectors can i ask a question smaller um and you hear epsilon q is is the or the line of jones well that's the only parameter i don't think you're among all the atom types in the head group somebody there in the room with that dress so we tried different questions protocols so we started off by using um bar kind of as a global metamodel and then um add some some kind of penalty function um if if the m bar arrows were too large but this doesn't work it didn't work so well for the reason that we've just seen that the rebending really is only um good in a very very close region and then the the errors will become small so it actually tricks you oh yeah yeah okay yeah yeah we had we had um like we had the m bar error we also had errors from different replica replicates and um but it still didn't work very well um and then what we ended up with is actually quite similar similar to force bounds so but we just get the local gradient from reweighting um we have the weights and the optimization that are chosen um by the user but also the uncertainty of the sensitivity plays a role in assigning these weights or in in modifying these weights in every iteration so that we only only um follow the directions where we kind of know um where we are very confident that we're going to see a a beneficial change in the parameters um can you guys hear us now no yeah this means both microphones if you accidentally tap it and i ask a question okay i think we're still waiting for them to be able to hear us them on but you can hear me but he says yeah that's about there and messaging them about it okay to take questions yeah so i think beno and michael shirts have questions can you hear us now hear me my do you hear my question can you turn them on here we're still working on them on or will you will you hear this testing and we have to go to your testing we must be muted in the room just a moment uh trying to get the audio turned up where's your here we go testing your audio testing output now we can hear you yeah okay testing okay now you got it great so maybe ben walk and ask this question i'll go after that can i talk okay uh you you talk about you talk about changing charges but i was not going to ask but now you you keep talking about charges what happened when you you change a charge do you keep you transfer like the opposite but equivalent charge uh increment to another group to keep like dipoles neutrals what do you do right yeah i should have mentioned that like we define charge groups um in the system to have um right to keep the overall charge neutral and to keep to keep the charge where do you put where do you put the the remainder the neighbors to the to the neighbors to the to the nearest neighbors uh in the charge group or to the to the other apartments in the charge group okay okay so when you say you change the charge of the carbonyl of the lipid what you're saying is that you changed the carbonyl the carbon and the oxygen at the same time right yeah we changed the oxygen directly and then the carbon changes indirectly exactly okay okay so follow up on that you showed certain molecules had larger smaller sensitivities did you look at the correlation between the charge magnitude and the change in the free energy what what i had observed you know sort of i don't think we ever published anything on this but there's a high correlation between how sensitive it was and how large the charge was larger charges changed by point zero you know by a delta um had a much bigger effect than smaller charges changed by the same delta okay this would be easy easy to pull up and check that yeah that would be good yeah okay thank you but and i'm very interested in where you have large can so in some of the messily papers there's the formula for the number of effective samples i'm catching up on questions i missed when the mute was happening uh and there essentially the way to think about it is if all the weights gets concentrated in one or two samples that leads to a low number of effective samples so there's there's a a way to convert that thinking in terms of density which weights are large into um in which are small into a number of effective samples it works quite well and we found that like a cutoff of 50 is where things start going haywire uh and then uh i'd be very interested on the ones like where epsilon has smooth behavior but the trends are quite different i'd be very interested in the future to figure out what's going on in those cases for sigma where the configuration space of us goes badly immediately that's sort of expected but it's it's actually interesting for epsilon that the trend is not quite right even though the uncertainty is low okay yeah all right yeah because there you you would expect um large confirmation or sufficient confirmation all from configuration over and up right and i think so i think rich messily did actually see the agreement was pretty good but that was for uh for homogeneous systems right yeah yeah so that's why'd be interested to see where the difference comes in if for some systems it works and some it didn't yeah yeah okay i've i've caught up with my questions i was asking over the last 10 minutes all right thank you um so yeah i mean actually the the when we did the free energy fitting uh then um we also didn't see the same like that there the sensitivities for epsilon were were top notch and here not like for the liquids not so much so that's that's definitely something to look into very interesting yeah um then let's let's so we did two iterations of this protocol um for all the atom types in the head group and um so this is the original c36 um with uh with uh just PME for charges and cut off monochones and uh so we haven't evaluated the the monolayer areas but we know from the surface tangents that they are really bad um and then if we use like but but all these these other things agree very well with experiment all all the other um areas um of the different biolayers um and then if we just use long range junks and we get um a large difference suddenly like this a five percent well more than five percent difference in the area polypids that decrease due to these original uh due to due to these attractive long range forces but after two iterations we're already back at the quality of the c36 force field and by um using when we look at at the monolayers the monolayers also um give very good match with experiment compress compressibility is a little worse um but but all the other biolayers that we tested so far look really good and this is after two iterations each iteration takes like three days we do a couple of simulations and um and it's kind of it was kind of surprising to us that by doing this this linear scheme you could get the proper or so such you can basically um capture the these non-linear um dependencies for a person observables and this would be really hard to do for a human without this kind of um without this kind of um tool i mean jeff did it but like i i wouldn't be able to do it i'm sure um so then when we look at the order parameters when we like those are the um the red dots are the experiments for the two chains and in the head group region and those are pretty good for the original c36 but then when we use ljpme it's all um it's all over the place and then again after two iterations we have the same we have basically the same accuracy as the original force field and in the head groups it's even better um right and i talked a little bit about having this regularization term where we basically have c36 as a target in the optimization um which is very moderately weighted but which prevents overfitting um and by that we waiting to see that the modifications on the final parameters are actually pretty minimal so charges change not more than 0.1 elementary charge units and the signals don't change more than 0.1 angstrom um and so far everything we've seen in terms of observables looks really good so far for the pc head group but we're going to do um the other other head groups too so i think we finally have a parameter set for the trial and before that works for long range on the domes and for the for monolayers and bilayers um and then i want to just spend maybe three or four minutes talking about where we're going with this so the first thing is that when we um when we when we have force fields that are editable like that's that's the dream right like we want a force field that we can that we can improve at any point in time but um as i said before one one um aspect of the charm force field that makes it popular is that it's been stable and that it's been good for many many decades or many years and uh and decades and um so how can we how can we get both like reliability for people to to trust in a force field but on the other hand keep it maintainable and keep um enable editing not like the tip-free water model that's like in there we're never going to get it out um and and i think one way to achieve this is continuous integration and i don't know how much you have thought about this but i think if if we develop force fields that we can edit in the future um then we have to think about how we how are we going to make people trust in the force fields um and one way to do this would just be a continuous integration framework where we have a repository where we save the final trajectories so that we can do reweighting on these trajectories easily and integrate new observables where we have the sensitivities saved so that if people want to add new simulations to parameterize against that we still have the sensitivities from the ultimulations to guide the optimization and prevent everything else or the previous optimization from being in vain and also have the fitting script all that online where we can basically ensure everybody can look at the results online everybody can use um our trajectories and our fitting procedure um but at the same time we make sure that the old or the the observables that we originally parameterized forward doesn't um don't screw up in the process i think this is an objective that we share for the open force field initiative in fact if we can assume some of the infrastructure burden for being able to do this for a wide variety of force fields through something like our property estimator which automates the computation of a lot of these observables maybe not the ones that you have yet but i think that would be fantastic to have that web framework that allows us to look at how are all of these different force fields doing on these kinds of physical quantities that are measured for and then um the the idea that you mentioned of being able to deploy the the trajectories as well would be great the trajectories can sometimes be very large so at least the large for a github repository but at least the ability to regenerate them and cache them locally if you want to do this exploration locally is something that we should support as well with this infrastructure yeah so i think we should definitely talk about how we can do enable that for the community yes okay and then finally just as um what just to say what we are going to do with this so um in the long term we want to talk we want to optimize also the jude force field the jude lipid force field is not in such a great state at the moment um in terms of the compressibilities for example and um we're going to try to use or we want to do to use this approach um for the polarizable charm force field as well so that we can have a longer term solution to improve the description of permeabilities and we also have this equivalent formulation of the jude model um multiple induced dipole model i should have put the citation here that's um Jing Huang published this two years ago which has basically um it's a translation of jude into multiple induced dipole framework but gives a little bit more flexibility in terms of the choice of parameters um and that's also something uh Bernie is very much uh interested in um so so these are going to be the things that we're going to work on um over the next couple of years um and your lunas is kind of play a big part in that and um with this i want to come to the conclusion so we first seen this automated approach to improve the description of the permeabilities but we can use multistage reweighting as kind of a global metamodal and then we have used similar methods to adapt c36 to linear zones pme um where we use reweighting as a local metamodel very much like force balance um where only in only two iterations we get to a force field that so far gives us everything that we want from from a from a little bit force field um but more things to test of course and that finally has consistent parameter settings for vitias and monolayers with this i um want to thank everybody who contributed to this word mostly a loon who has done a big part of the uh lipid refitting uh project um and thank you all for your attention take questions does anyone who questions like in r i think i have a question yeah please do just just to uh to really understand what happened so when you run the old charm force field with the sort of 12 ish instrument cut off the area per lipid let's say for a dppc was something like 63 inches square i think 62.9 that's what you have in the table i think and then when you turn on this uh pme linear zones it dropped to about 58 that's right and then by adjusting linear zones and partial charges for basically any atone or just a few atones i don't that's so sure i mean how many parameters were changed you basically go back up to about 63 angstrom that is right yeah we we train on every atom in the head group um and the reason for this is that most of these parameters are very heuristic like the partial charges um are not like the the derivation of the partial charges in the head group is not very systematic it's it's a little bit of trial and error um for many of these charge charges and um one other thought a line of thought was that when we distribute the burden of of capturing this um overall parameters in the head group then we can get away with very small changes and thereby hopefully also ensure compatibility with the protein force field for example and i think that like from what we see in the parameters that they that they change so little the final parameters we can hope that we still get um good interaction between the lipid and the force in the protein force field for example let's see so you changed also the carbons between the uh glycerol backbone and the phosphate like there are a couple of aliphatic carbon type you don't change just the phosphate and the uh you know the very polar group right you change like all the carbons in between yeah we change the link atoms between the phosphate and the and the choline for example yeah okay but you don't change the dihedral or you adapt the dihedral yes we do we do adapt the dihedral um so we have we tried out different approaches for this what we're going for now is that a lot of the dihedral distributions we fit to the c36 distributions because we know that they are pretty good um and that's a very doable approach and there's also in the spirit of not deviating too much from the original force field um so we're readjusting the dihedral's to give the same distributions as c36 if possible i mean distribution in the simulation right in the in the bio simulation so okay okay so you calculate like the distribution of the dihedral in the c36 with 12 inch from cutoff and then you try to maintain those dihedral distribution exactly except and except in the cases where we definitely know a better answer from the sd or from qm and we have a trustworthy or more trustworthy reference um than c36 except we like but but for most of the dihedral's we we fit them to c36 i see i see it may be a minute it's very interesting the um it may be that the the uh the tweaking of partial charges would have to be looked at in a maybe slightly different light in the case of the druid force field because all the charges have been fitted from ab initio to begin with yeah small change probably are still acceptable but they they're not as heuristic as charm 36 charges right so so it would be like one first step for the druid force field would be to look at the reweighting for both the polarizabilites and for the charges and to see like which is more which of these two is more sensitive and what do we um like if the one is heuristic then the other is not right so i believe in in the in the lipid force field you you did the uh chart the partial charges from ab initio that but then you had heuristic information for the polarizabilites right well the polarizabilites came from ab initio too but they they are a little bit more uh they're adjusted more empirically to some extent too at the end of the day like you know to get dietary constant right and things like that so right the the the charges are a little bit more restrained by the ab initio the polarizabilites a little bit less they're also extracted from ab initio but if we can get away with with changes in the same order of magnitude like for example point one electric start chart elementary charge unit that's like um probably also acceptable for something that is fit to ab initio yeah just just to finish here it's i think that we highly suspect that the uh incompressibility of the bilayer with the druid force field comes from the polar red polar red interaction and uh that probably would have required doing some sort of basmatic pressure measurement or something like that to really quantify how much the like for example the choline binds to the phosphate and things like that this was this was a bit unanticipated but this kind of problem between charged species has has been observed more regularly like between arginines and aspartic acid like a lot of uh the very charged stuff has required a second look afterward mm-hmm yeah no okay okay thank you but we're we're definitely going to be in touch about about this effort at least any other questions well there are no questions here in the room i now we can hear you on the room i apologize for that i checked you like if you can hear us but not if you can hear you so um sorry about that but if there are no questions i think we can finish with you thank you so much and dress and thank you