question regarding epik output files

[DivyaMatta15]

Just a quick question. I was trying to understand the output files, its bit challenging for me. How does the numbers tells the connectivity of atoms(how to know which atom connects to which atom) for example, in histidine what 1 1 2 1 means ? Also what does DICT means?

Thanks

[jchodera]

Which file(s) specifically are you talking about?

[jchodera]

The ACE-HIS-NME example is here. You should be focusing on two different output files:

• The mol2 file is a standard Tripos mol2 format documented here that contains atom names, partial charges, and bonds. There are many readers for Tripos mol2 files, and the example script we had written last year has an example of mol2 file reading using the OpenEye toolkit.
• The state penalties file lists the relative free energy of each of the protomer states (in kcal/mol) at pH 7.4. These should be used to compute the effective penalties or pKas you use for a simulation only at pH 7.4.

[DivyaMatta15]

I was looking at the histidine output mol2 file. The document is helpful. I am looking into it.

Thank you!

[newbooks]

I am looking at ACE-HIS-HIS-NME.mol2. It contains one charge set so far. Do we have charges at other ionization states?

Also where is supposed be the chain ID column, if you anyone happens to know? I will read the documentation as well. So no worry on this.

[newbooks]

It looks like mol2 ATOM records do not include chain ID. Correct me if I am wrong.

[jchodera]

I am looking at ACE-HIS-HIS-NME.mol2

You are looking at the input files. Go here to find the output files:

https://github.com/choderalab/mcce-charges/tree/master/his-test-peptides/output/HIS-HIS

You should focus on the output mol2 file here:

https://github.com/choderalab/mcce-charges/blob/master/his-test-peptides/output/HIS-HIS/HIS-HIS-epik.mol2

MOL2 files do not contain chain IDs.

[mrgunner]

For the diHis there should be 9 conformers. 4 with 0 charge, 4 with +1 and 1 with +2 (Divya can enumerate). We are hoping to see epic can provide these so formers. If it won't I'd like to understand what it is doing.

Just add the chain designator (I think MCCE will do this automatically)

On Apr 18, 2016, at 3:42 PM, Junjun Mao notifications@github.com wrote:

I am looking at ACE-HIS-HIS-NME.mol2. It contains one charge set so far. Do we have charges at other ionization states?

Also where is supposed be the chain ID column, if you anyone happens to know? I will read the documentation as well. So no worry on this.

— You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub

[SalahBioPhysics]

Looking at the epik output I can see 9 sections (perhaps that's the 9 conformers) starting with atTRIPOSATOM. However, as Dr. Gunner mentioned we need 4 conf with 0 charge, 4 conf with +1 and 1 conf with +2 charge. Where in the epik output we can find these different states?

Right now, we will use the epik output to generate the tpl files needed to run mcce. After we do the di-his, and test it, we can run the KINASE INHIBITORS.

[jchodera]

Please look at the mol2 output files I prepared for you along with the state penalties txt file. The mol2 file has all of the protomers in a standard mol2 format, along with charges. The total charge on each molecule should add up to the integral charge of each protomer.

You should not need to look at the Epik output files. We have already extracted that information for you.

[newbooks]

Is this file the output? https://github.com/choderalab/mcce-charges/blob/master/his-test-peptides/output/HIS-HIS/HIS-HIS-epik.mol2

I browsed it and it seems to me the HIS2 and HIS3 in entry 2 and 5 essentially have identical atoms, coordinates and charges. They differ by extra small atom H46 and H47, which have 0 charges (9th column) and they don't make difference in continuum model.

Any opinions?

Junjun

On Wed, Apr 20, 2016 at 7:43 PM, John Chodera notifications@github.com wrote:

Please look at the mol2 output files I prepared for you along with the state penalties txt file. The mol2 file has all of the protomers in a standard mol2 format, along with charges. The total charge on each molecule should add up to the integral charge of each protomer.

You should not need to look at the Epik output files. We have already extracted that information for you.

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-212653834

[jchodera]

You want the -charged.mol2 suffix: https://github.com/choderalab/mcce-charges/blob/master/his-test-peptides/output/HIS-HIS/HIS-HIS-epik-charged.mol2

[SalahBioPhysics]

YES, thanks a lot, that was useful

[jchodera]

No problem. My fault for not creating a README.md indexing the files. I'll be sure to make one!

[newbooks]

The atom names in "charged" output do not match another mol2 file. Is it safe to assume an one on one name mapping? Or I should match the names by their coordinates to be correct?

Salah, "charged" mol2 file do not contain residue information. I guess you will need to load both mol2 files to get things MCCE need.

Junjun

On Wed, Apr 20, 2016 at 9:50 PM, John Chodera notifications@github.com wrote:

No problem. My fault for not creating a README.md indexing the files. I'll be sure to make one!

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-212692156

[jchodera]

The atom names in "charged" output do not match another mol2 file.

Which "another mol2 file" are you referring to? Is this important?

The multiple protomers within the -epik-charged.mol2 file all have the same heavy atom names, though their proton names may differ. I didn't think this would be problematic.

I am under the impression that you just need to create a discrete protomer library for a single conformer. When we create a TPL file for the small molecule ligand, it will be in a single conformation fixed in space.

[jchodera]

If you can give us more insight into the problem you are having, we may be able to provide a bit more help here. I just don't understand what the issue you're having is.

[newbooks]

I am referring https://github.com/choderalab/mcce-charges/blob/master/his-test-peptides/output/HIS-HIS/HIS-HIS-epik.mol2

and https://github.com/choderalab/mcce-charges/blob/master/his-test-peptides/output/HIS-HIS/HIS-HIS-epik-charged.mol2

The names in mol2 is what we want, but as you said we should use the charge in charged.mol2. We are mapping the names.

In mol2 file: 1 HH31 2.2108 -0.3185 0.7166 H 1 ACE 0.0640 2 CH3 1.1512 -0.0260 0.7573 C.3 1 ACE -0.1731 3 HH32 0.6861 -0.4570 1.6561 H 1 ACE 0.0640 4 HH33 0.6331 -0.3991 -0.1385 H 1 ACE 0.0640 5 C 1.0475 1.4623 0.8084 C.2 1 ACE 0.6562 6 O -0.0576 1.9895 0.8580 O.2 1 ACE -0.5801 7 N 2.2668 2.1122 0.7917 N.am 2 HIS2 -0.5529 BACKBONE|DICT|DIRECT 8 H 3.1660 1.4804 0.7451 H 2 HIS2 0.3087 DICT 9 CA 2.3956 3.5471 0.8321 C.3 2 HIS2 0.0354 BACKBONE|DICT|DIRECT 10 HA 1.5707 3.9122 1.4616 H 2 HIS2 0.0955 DICT 11 CB 3.7169 3.9924 1.4694 C.3 2 HIS2 -0.1009 DICT 12 HB2 3.9486 3.3418 2.3256 H 2 HIS2 0.0696 DICT 13 HB3 4.5247 3.9218 0.7261 H 2 HIS2 0.0696 DICT 14 CG 3.6731 5.4156 1.9697 C.2 2 HIS2 0.2849 DICT

In charged.mol2: 1 H1 2.2108 -0.3185 0.7166 H 1 <0> 0.0661 2 C1 1.1512 -0.0260 0.7573 C.3 1 <0> -0.1720 3 H2 0.6861 -0.4570 1.6561 H 1 <0> 0.0661 4 H3 0.6331 -0.3991 -0.1385 H 1 <0> 0.0661 5 C2 1.0475 1.4623 0.8084 C.2 1 <0> 0.6798 6 O1 -0.0576 1.9895 0.8580 O.2 1 <0> -0.5994 7 N1 2.2668 2.1122 0.7917 N.am 1 <0> -0.5767 8 H4 3.1660 1.4804 0.7451 H 1 <0> 0.3237 9 C3 2.3956 3.5471 0.8321 C.3 1 <0> 0.0503 10 H5 1.5707 3.9122 1.4616 H 1 <0> 0.0953 11 C4 3.7169 3.9924 1.4694 C.3 1 <0> -0.0206 12 H6 3.9486 3.3418 2.3256 H 1 <0> 0.0704 13 H7 4.5247 3.9218 0.7261 H 1 <0> 0.0704 14 C5 3.6731 5.4156 1.9697 C.2 1 <0> -0.2269 15 N2 4.0317 6.4557 1.1492 N.pl3 1 <0> -0.3071 16 C6 3.8805 7.5324 1.8880 C.2 1 <0> 0.4015 17 H8 4.0690 8.5501 1.5765 H 1 <0> 0.0622 18 N3 3.4415 7.2272 3.1463 N.2 1 <0> -0.6676 19 C7 3.3031 5.8694 3.2137 C.2 1 <0> 0.3052 20 H9 2.9660 5.3709 4.1112 H 1 <0> 0.0378

On Thu, Apr 21, 2016 at 3:31 PM, John Chodera notifications@github.com wrote:

If you can give us more insight into the problem you are having, we may be able to provide a bit more help here. I just don't understand what the issue you're having is.

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213081603

[jchodera]

Why are you mapping the names between the two? Why do you need both files?

We can fix the naming inconsistency if that is important, but I don't understand why you are trying to use both files.

[newbooks]

Yes, as you can see with HIS-HIS-epik-charged.mol2 file alone, we don't have residue information and atoms seem to be named sequentially (I am afraid this will cause trouble if an input file has atoms in a different order).

We already have a script to fix it though. A script will extract charges from charged.mol2 and place that in mol2 file.

Junjun

On Thu, Apr 21, 2016 at 9:47 PM, John Chodera notifications@github.com wrote:

Why are you mapping the names between the two? Why do you need both files?

We can fix the naming inconsistency if that is important, but I don't understand why you are trying to use both files.

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213207231

[SalahBioPhysics]

At that note I want to mention that we don't have epik-charged.mol2 files for the epik_inhibitors output.

[newbooks]

Let's work on standardizing the parameterization of HIS-HIS. Once a procedure is laid out, getting more inhibitors should be easy.

I am away next week. Will keep looking after I am back on Friday.

Junjun

On Fri, Apr 22, 2016 at 2:18 PM, Salah Salah notifications@github.com wrote:

At that note I want to mention that we don't have epik-charged.mol2 files for the epik_inhibitors output https://github.com/choderalab/mcce-charges/tree/master/epik_inhibitors/output .

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213539638

[andrrizzi]

At that note I want to mention that we don't have epik-charged.mol2 files for the epik_inhibitors

Yeah, sorry that's my fault. I can easily create them after we have worked out the HIS-HIS case. Or sooner in case you need them ASAP.

[SalahBioPhysics]

No, don't need them now, I'm just giving you heads up :)

[newbooks]

Hi all,

We have a bit of methodology incompatibility issue here.

It looks to me epik calculates the charge for the whole molecule, in our case, it's all 4 residues in ACE-HIS-HIS-NME. Therefore, an individual residue doesn't necessarily have an integer charge. For example 6 atoms in the first residue ACE gives 0.0661-0.1720+0.0661+0.0661+0.6798-0.5994 = 0.1067

MCCE on the other hand treats a residue as a separable unit, which must have an integer charge.

Any ideas?

Junjun

On Fri, Apr 22, 2016 at 2:18 PM, Salah Salah notifications@github.com wrote:

At that note I want to mention that we don't have epik-charged.mol2 files for the epik_inhibitors output https://github.com/choderalab/mcce-charges/tree/master/epik_inhibitors/output .

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213539638

[jchodera]

JunJun, we decided long ago that for now, we will treat small molecules as A SINGLE RESIDUE. We can explore ways to automatically partition molecules into independent residues in the future, but this will be a major research project.

All you need to do is modify the TPL generation script to convert the whole molecule into a single residue with a static conformation and multiple protomers. You will need to calibrate the model pKas or reference energies to reproduce the Epik-computed populations (specified as state relative energies in kcal/mol) when the residue is simulated at the pH used to generate these populations (here, pH 7.4).

[newbooks]

John, yes, as far as I know, we will do one as single molecule and another as separate residues as control, using a well known system HIS-HIS.

An quick fix is normalize the atom charges residue by residue. That should be close enough for now. Can other people weigh in?

Junjun

On Fri, Apr 22, 2016 at 5:53 PM, John Chodera notifications@github.com wrote:

JunJun, we decided long ago that for now, we will treat small molecules as A SINGLE RESIDUE. We can explore ways to automatically partition molecules into independent residues in the future, but this will be a major research project.

All you need to do is modify the TPL generation script to convert the whole molecule into a single residue with a static conformation and multiple protomers. You will need to calibrate the model pKas or reference energies to reproduce the Epik-computed populations (specified as state relative energies in kcal/mol) when the residue is simulated at the pH used to generate these populations (here, pH 7.4).

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213600969

[mrgunner]

Junjun It's great to have you on board. As john says we will do the inhibitors as multiple conformers not breaking them up as we usually do (like we do with the heme to divide and conquer). The aim with the diHis is to make sure we get the right (I think 9) conformers with epic ((1w/crg 2) and 4 each with crg 1 and 0 ) and that the energies of the states are close to what we would get with mcce ( here there are 2 residues each with 3 conformers -giving the same 9 micro states). The individual charges for the residues in the epic conformers will not likely be integer, thou the sum over the whole diHis conformer will be integer.

My view has been if we can translate the epic output for the diHis we will be able to handle (and trust) the inhibitors, where we don't know exactly what to expect.

It sounds like progress is being made. M

On Apr 22, 2016, at 4:14 PM, wJunjun Mao notifications@github.com wrote:

Hi all,

We have a bit of methodology incompatibility issue here.

It looks to me epik calculates the charge for the whole molecule, in our case, it's all 4 residues in ACE-HIS-HIS-NME. Therefore, an individual residue doesn't necessarily have an integer charge. For example 6 atoms in the first residue ACE gives 0.0661-0.1720+0.0661+0.0661+0.6798-0.5994 = 0.1067

MCCE on the other hand treats a residue as a separable unit, which must have an integer charge.

Any ideas?

Junjun

On Fri, Apr 22, 2016 at 2:18 PM, Salah Salah notifications@github.com wrote:

At that note I want to mention that we don't have epik-charged.mol2 files for the epik_inhibitors output https://github.com/choderalab/mcce-charges/tree/master/epik_inhibitors/output .

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213539638

— You are receiving this because you commented. Reply to this email directly or view it on GitHub

[jchodera]

I wouldn't do that. That is not what we plan to do for actual kinase inhibitor molecules, so it would be useless as an exercise to test how well your pipeline will work for kinase inhibitors.

[newbooks]

It all makes sense if we don't do control. For a new cofactor, what we need are 1) atom charges 2) solution pKas if the cofactor has protonation states

The epik output gives out 9 states and corresponding charge sets. Presumably solution pKas can be derived from the state penalties.

I need to understand a bit more about how epik state penalties are calculated. To be more specific, at what pH and what medium dielectric constant are the state penalties in epik calculated? Or, how do I convert state penalties to state population at various pHs as in the graph we saw at John's office.

Thanks,

Junjun

[jchodera]

How does MCCE compute protomer acceptance probabilities at a given pH?

Suppose we have a single residue in solution. If the acceptance criterion gives a probability

P_i = exp[-beta(U_i + ln(10) (pH - pKa_i))] / Z

then we need to determine the pKa_i artificially to give us

P_i = exp[-beta (state_penalty_i + c) ]

where c is some constant.

Is there another mechanism besides assigning pKas to individual protomers that can be used to adjust their equilibrium populations?

[newbooks]

How does MCCE compute protomer acceptance probabilities at a given pH? MCCE uses Monte Carlo sampling. It determines the acceptance of a new state by Boltzmann distribution equation and the energy difference between the new state and the current state. The state energy includes self energy and pairwise interactions. pH effect is part of self energy calculated similarly as your equation.

Is there another mechanism besides assigning pKas to individual protomers that can be used to adjust their equilibrium populations? Yes, we reserved EXTRA term, whose original purpose was to include energy contribution from un-accounted sources, in the form of self energy of each conformer.

Due to epik state penalties acting as a "total" state energy term, they may not be directly incorporated in MCCE as a constant self energy term. MCCE divides self energies into several types and some of them varies with dielectric boundary.

Is it easy for your group to write a function to convert epik state penalties to state population at a given pH? We will parameterize the inhibitor to reproduce that titration curve of single molecule in solution.

Junjun

On Fri, Apr 22, 2016 at 10:13 PM, John Chodera notifications@github.com wrote:

How does MCCE compute protomer acceptance probabilities at a given pH?

Suppose we have a single residue in solution. If the acceptance criterion gives a probability

P_i = exp[-beta(U_i + ln(10) (pH - pKa_i))] / Z

then we need to determine the pKa_i artificially to give us

P_i = exp[-beta (state_penalty_i + c) ]

where c is some constant.

Is there another mechanism besides assigning pKas to individual protomers that can be used to adjust their equilibrium populations?

— You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/choderalab/mcce-charges/issues/17#issuecomment-213651367

[jchodera]

MCCE uses Monte Carlo sampling. It determines the acceptance of a new state by Boltzmann distribution equation and the energy difference between the new state and the current state. The state energy includes self energy and pairwise interactions. pH effect is part of self energy calculated similarly as your equation.

Yes, but we need to know exactly what stationary distribution is sampled in terms of the parameters you can assign each protomer in MCCE. For example, how do the pH, pKa for each tautomer, and EXTRA terms enter into either the stationary distribution or acceptance probability?

Is there another mechanism besides assigning pKas to individual protomers that can be used to adjust their equilibrium populations? Yes, we reserved EXTRA term, whose original purpose was to include energy contribution from un-accounted sources, in the form of self energy of each conformer.

This might be what we want to use, then. Can you point me to the place in the manual that describes how this enters into the stationary probability?

Due to epik state penalties acting as a "total" state energy term, they may not be directly incorporated in MCCE as a constant self energy term. MCCE divides self energies into several types and some of them varies with dielectric boundary.

Is it easy for your group to write a function to convert epik state penalties to state population at a given pH? We will parameterize the inhibitor to reproduce that titration curve of single molecule in solution.

It is important to remember that we are only trying to work at a SINGLE pH. The epik state penalties are generated assuming pH 7.4 is the only pH we ever care about. We can compute whatever you need, but we need to know what it is you need first.

If you can explain how the EXTRA term is used along with the potential energy of each conformer, that would help us figure out what you need. Also, is the internal energy of interactions within the residue included in the energy, or just interactions between residues?

[newbooks]

I try to summarize the procedure in a few sentences:

1. Protein is broken into residues. Each residue has a backbone and side chain conformers which differ in charge state or side chain position.
2. Each conformer has energy G = Gbackbone_ele + Gbackbone_lj + Gtorsion + 2.3_m_kbT*(pH - pK0) + Gextra + dGrxn + Sum(Gpairwise_ele + Gpairwise_lj)
3. A microstate is a selection of one conformer from each residue, and the state energy is the sum of above individual conformer energy.
4. Monte Carlo sampling is then performed to get acceptance of microstate based on microstate energy. Microstate is traced but not recorded. The statistics of conformers in accepted microstate is recorded and used to calculate protonation.

Gbackbone_ele: Electrostatic interaction to backbone Gbackbone_lj: Lennard Jones interaction to backbone Gtorsion: Torsion energy 2.3_m_kbT*(pH - pK0): pH effect, m is 1 for base or -1 for acid. pK0 is solution pKa of the residue and kbT is a constant. Gextra: Extra energy term to scoop energy from all un-accounted sources, usually small for amino acids, as we have good understanding to these. dGrxn: Change of reaction field energy of this conformer in water and in protein. Sum(Gpairwise_i): Pairwise interaction to other selected conformers in other residues.

EXTRA term (MCCE coders: this might be named as OFFSET in the code) and pK0 play the same role and are complimentary. if pK0 is unknown, it can be accounted in EXTRA.

We can parameterize the term (-2.3_m_kbT*pK0 + Gextra) if we know the state population at a single pH. For a single molecule in water, dGrxn + Sum(Gpairwise_ele + Gpairwise_lj) is 0.

Junjun

[jchodera]

Thanks for the detailed summary!

I presume that we can specify Gextra for each conformer. If that's not correct, let me know and I can propose a different idea.

Here's what I would suggest:

• If it is easy to do so, set m = 0 so that there is no contribution from the pH effect term 2.3*m*kbT*(pH - pK0). Since the Epik state penalties change in a complex way with pH, there is simply no way for us to transform those energies into site pKas that would make sense here. We will put the state penalties into another term. if you can't set m = 0, set pH = pK0 so the contribution is still zero.
• Set dGrxn = 0. We won't use this either.
• As long as all protomers of the new residue have the same torsion energy Gtorsion, it doesn't matter what this is since we are only dealing with one conformation. If the Gtorsion energies for the different protomers are different because they have different numbers of torsions, we will have to subtract off those contributions from the Epik state penalties.
• If there is no intraresidue energy contribution to Gbackbone_ele, Gbackbone_lj, Gpairwise_ele, and Gpairwise_lj, then we put the whole Epik state penalty for each protomer into the Gextra term for that protomer`.

If we do this, then the total energy for just the residue in solution should be

G = Gbackbone_ele + Gbackbone_lj + Gtorsion +
2.3*m*kbT*(pH - pK0) + Gextra + dGrxn + Sum(Gpairwise_ele + Gpairwise_lj)
= 0 + 0 + 0 +
0 + Gextra + 0 + Sum(0 + 0)
= epik_state_penalty

This would mean that we achieve the correct distribution of protomers for the residue alone in solution for pH 7.4---which is what Epik reports. When we put this residue into a protein context, the electrostatic and LJ interactions would perturb the state population.

Let me know if I've misunderstood this and you can't set Gextra or pK0 separately for each protomer.

[mrgunner]

I was at my family yesterday but the proposed treatment will work. Turn off internal energies and put the epic energy into extra. (And set nH to 0). We can only do the calc at one ph.

Not now but in the future it oils be interesting to see epic calc as a f(ph). But afterwards.

M

On Apr 23, 2016, at 8:52 AM, John Chodera notifications@github.com wrote:

Thanks for the detailed summary!

I presume that we can specify Gextra for each conformer. If that's not correct, let me know and I can propose a different idea.

Here's what I would suggest:

If it is easy to do so, set m = 0 so that there is no contribution from the pH effect term 2.3_m_kbT*(pH - pK0). Since the Epik state penalties change in a complex way with pH, there is simply no way for us to transform those energies into site pKas that would make sense here. We will put the state penalties into another term. if you can't set m = 0, set pH = pK0 so the contribution is still zero. Set dGrxn = 0. We won't use this either. As long as all protomers of the new residue have the same torsion energy Gtorsion, it doesn't matter what this is since we are only dealing with one conformation. If the Gtorsion energies for the different protomers are different because they have different numbers of torsions, we will have to subtract off those contributions from the Epik state penalties. If there is no intraresidue energy contribution to Gbackbone_ele, Gbackbone_lj, Gpairwise_ele, and Gpairwise_lj, then we put the whole Epik state penalty for each protomer into the Gextra term for that protomer`. If we do this, then the total energy for just the residue in solution should be

G = Gbackbone_ele + Gbackbone_lj + Gtorsion + 2.3_m_kbT*(pH - pK0) + Gextra + dGrxn + Sum(Gpairwise_ele + Gpairwise_lj) = 0 + 0 + 0 + 0 + Gextra + 0 + Sum(0 + 0) = epik_state_penalty This would mean that we achieve the correct distribution of protomers for the residue alone in solution for pH 7.4---which is what Epik reports. When we put this residue into a protein context, the electrostatic and LJ interactions would perturb the state population.

Let me know if I've misunderstood this and you can't set Gextra or pK0 separately for each protomer.

— You are receiving this because you commented. Reply to this email directly or view it on GitHub

[jchodera]

Not now but in the future it oils be interesting to see epic calc as a f(ph). But afterwards.

This is easy to automate, so it should be no problem if we can get the proposed workflow going. But for right now, we just need to see what happens at pH 7.4.

[SalahBioPhysics]

I'm looking into the DB8.tpl.pH-7 file that was created with tpl_maker_am1bcc.py. Is it possible to add the output files for standalone epik? specifically the -epik.mol2 file and (if possible) the epik-charged.mol2

[SalahBioPhysics]

Can I use the output from bosutinib?

[jchodera]

The bosutinib output you link is the equivalent output with updated settings for Epik. If you want to try to convert that directly, that is what you would use. That particular output doesn't correspond exactly to the older DB8.tpl.pH-7 file because we don't feed Epik the same settings, though.

If you let me know more about what you're having trouble with or would like to do, I might be able to be of more help. I think the best strategy is for you to modify the tpl_maker_am1bcc.py to use the *-epik-charges.mol2 files we generated for you and produce output in your current TPL format.

[jchodera]

Actually, it looks like that directory doesn't have -epik-charges.mol2 files because we haven't modified the script that generates these mol2 files to correctly assign charges to the Epik-enumerated protomers. We already made the change to the corresponding script in his-test-peptides, so let me make the change to the kinase inhibitors script and regenerate these kinase inhibitor charged mol2 files.

[jchodera]

@SalahBioPhysics : I've now updated the kinase inhibitors for you, adding the charged protomers and state penalty text files in https://github.com/choderalab/mcce-charges/pull/18.

See the updated bosutinib directory, with these files in particular:

• Bosutinib-epik-charged.mol2 - use this file to get the different protomers and their charges
• Bosutinib-state-penalties.out - use this file to get the relative free energies of each corresponding protomer (in kcal/mol) for the inhibitor in solvent

All the corresponding files for other inhibitors are also available in the output directory.

[jchodera]

@bas-rustenburg : You should be able to use these input files for the constant-pH simulations too.

[SalahBioPhysics]

Thanks, yes that's what I'm doing right now, modifying the tpl_maker_am1bcc.py to use the *-epik-charges.mol2 files.

[SalahBioPhysics]

Does epik uses a PDB file as an input? If yes, is it possible to give us a list of the PDB IDs used or upload them into here?

[jchodera]

I'm not sure what you mean. Epik parameterizes just the small molecule, and the small-molecule structures were generated from canonical isomeric SMILES strings. There were no PDBs involved here.

We could generate PDB files of the small molecules if this would be helpful for you. Just let me know---it's easy for us to do this!

[jchodera]

@SalahBioPhysics: Did I understand correctly that you want us to generate PDB files for the small molecules? If so, it's easy for us to generate them. Just let me know.

[SalahBioPhysics]

Yes, that would be good as we can use the connectivity section at the bottom of the pdb.

[jchodera]

OK, I'll generate these. Note that the mol2 files also have connectivity information you can use!

[jchodera]

Just to clarify: Do you want the PDB files to contain just one protomer, or all protomers (as a multi-MODEL PDB file)?

[SalahBioPhysics]

Thanks. One protomer. I'm modifying Denise tpl_maker.py program to read both .mol2 file and pdb. I could make it only read .mol2 files, but then I have to reinvent the wheel for the connectivity part as Denise program already does the connectivity but for an pdb input.