calculate thermodynamics for universal reactions / metabolites

[pstjohn]

This is obviously a big task, but I thought it was worth opening up an issue to discuss. It would be cool if there was a built in method to calculate and/or aggregate thermodynamic info on bigg reactions and metabolites, similar to equilibrator for KEGG and von Bertalanffy for the cobra MATLAB toolbox.

Metacyc also has a gibbs free energies for their database of reactions and metabolites, but I'm not sure the metanetx cross-references are enough to get great coverage of the bigg reaction database. component-contribution might be a good place to start, but it has many dependencies (including the academic-only Marvin for pKa calculations), and will take some work to include non-KEGG metabolites.

[jlerman44]

Does the von Bertalanffy for the cobra MATLAB toolbox package extend beyond KEGG and the component-contribution code of Elad/Ronan?

[jlerman44]

Make sure not to combine gibbs free energies calculated using different underlying models/assumptions.

[cdanielmachado]

Converting reactions from BiGG nomenclature to KEGG nomenclature and calculating with eQuilibrator is an option. But I've found a few limitations while doing this:

• BiGG compounds without associated KEGG id
• BiGG compounds associated with multiple KEGG ids
• Unbalanced BiGG reactions

I believe the first problem is by far the main obstacle. The solution for this, like Peter mentions, is to calculate formation energies for non-KEGG compounds using the component contribution code (which requires installing all the dependencies like Marvin and so on).

If someone can run this (just has to be done once) for all BiGG compounds, and add it to the database, then the reaction thermodynamics can be trivially computed for all reactions.

[pstjohn]

So I believe the matlab toolbox extends beyond KEGG, because they want you to have the cxcalc Marvin program installed.

@cdanielmachado, apparently computing free energies for reactions and compounds are somewhat different (and delta G's of whole reactions are much more accurate). Indeed, there are a number of compounds in BIGG that are not in KEGG. So we'd need someone with a Marvin academic license to run them for all the BIGG reactions. To do that, we'll need Inchi's or SMILES for all the bigg compounds, which is also not an easy task. (I can't believe there isn't an open source pKa calculator... that seems like such a basic QSPR, and there's no shortage of commercial options)

There's slightly better coverage of the BIGG metabolites / reactions in MetaNetX to MetaCyc, so that could also be a source of structural information.

[phantomas1234]

https://github.com/eladnoor/component-contribution also uses chemaxon for pKa predictions and also goes beyond KEGG. We're thinking about replacing the proprietary chemaxon dependency with rdkit.

[cdanielmachado]

Good points. Also, I forgot to mention one extra issue. For reactions involving compounds in more than one compartment (which is a large fraction of all reactions in BiGG), one must take into account the energy associated with moving compounds across compartments, which will depend on compound charge and membrane potential.

Without that, you easily generate unlimited energy from transport loops.

So, I take back my previous claim. Compound formation energies will not be sufficient to calculate reaction thermodynamics.

[pstjohn]

@phantomas1234, is there a function in rdkit that predicts a pKa? I wasn't aware of one, but that would be awesome. Or would you fit your own QSPR? You could also swap out the openbabel dependency for rdkit in either case (just to keep the requirements smaller).

Also, while component-contribution can extend beyond KEGG, the integration is pretty tight. Check out https://groups.google.com/forum/#!topic/equilibrator-users/Cli5xfQcVM0 for some additional considerations on extending to molecules outside of KEGG.

Let me know if you do try to move forward with this, I'm happy to contribute!

EDIT: Also just to make it explicit, thermo calculations in the MATLAB COBRA toolbox and the component-contribution repo are closely related projects and seem to share a lot of the same code.

[draeger]

Once we have such values, it would also be good to directly store them in the SBML files for download. At the moment, there is no field where this could be saved. The FBC package would certainly be the right place to store such information, which will require some discussion with the package working group.

[cdanielmachado]

@draeger: While I do like that option, I wonder how feasible it is on the long run to keep extending the fbc package with new reaction/metabolite attributes that people come up with.

For instance, some people may want to have ∆Gº values, others might also like to have the ∆Gº std values. Some people might say they want EC numbers to be reaction attributes as well.

Maybe it would be easier to just agree on storing optional attributes as key, value pairs in the annotations field? (At least that's how I'm doing it for the moment.)

[draeger]

Well, EC numbers can already be annotated using MIRIAM, so that's not a problem. I think, what we need here is some kind of an (at least semi) controlled key-value pair list, where people could request new keys as they become necessary and then attach a value to it.

For storing the standard reaction Gibbs energy, a parameter object could be used already. The only problem is that there is no direct link that would specify the relation between the parameter and the reaction it refers to. All we would need is a key-value pair with something like STANDARD_REACTION_GIBBS_ENERGY and the id of a parameter object attached to the reaction.

A trick could be to define the Keq value of each reaction using Keq = e^(-∆Gº/RT) and reuse this Keq in the kinetic law by defining the mass-action rate law as v = k+ (reactants - products/Keq). This would also give a link, which is much more difficult to unravel though.

[jlerman44]

@phantomas1234 Is the rdkit option ready for primetime?

[phantomas1234]

@jlerman44 turned out that rdkit was a dead end for pKa predictions and that there are only other commercial solutions that could do the job. Miguel knows more about this. Regarding the FBC package, I think it is a shame that thermodynamic FBA is not used more widely. Hope we can change that and it would be good to find a solution for SBML to store the necessary information.

[draeger]

Well, it depends on what you need to store in SBML. You can always add thermodynamic constants in form of local or global parameters, optionally with constant flag set to true. Using proper annotation, it should be clear what is meant. If you think, something more specific needs to be defined, please contact the FBC mailing list or sbml-discuss@googlegroups.com or at least e-mail to the SBML editors: sbml-editors@googlegroups.com.

[phantomas1234]

Hmmm, I don't know. What is the point of defining a normal SBML parameter if there is no corresponding equation? I guess it is possible to put the information in somehow but different tools might encode the information differently and then interoperability doesn't work. I was planning to get in touch with the FBC folks at some point.

[draeger]

Yes, this makes sense. This kind of information could either go into the annotation section or stored within a specific field where it is relevant. The FBC group is the right place to request such changes.

[djinnome]

Just so folks are aware, thermodynamic calculations were performed in the Feist et al 2007 iAF1260 publication.

[wbryant]

For what it's worth I used my own hacked-together version of eQuilibrator to predict reversibility of all KEGG reactions for use in fastGapFill. When I checked through my results manually against the reversibility in KEGG the eQuilibrator results were inconsistent in a significant fraction of cases. I don't know if this is useful/relevant, but I thought I'd mention it.

Also, we don't have an institutional subscription to KEGG where I am and I failed to determine an automatic way to scrape reversibility from the website - which is why I used eQuilibrator in the first place.

[cdanielmachado]

@wbryant how do you get reversibility information from KEGG?

[draeger]

I just posted the link to this discussion to the FBC package development group so that ways can be identified to store thermodynamic properties in a standardized way in SBML documents.

[rmtfleming]

It would be good if these could be followed: Recommendations for Nomenclature and Tables in Biochemical Thermodynamics DOI: 10.1111/j.1432-1033.1996.0001h.x

[wbryant]

@cdanielmachado

1. Find the relevant reaction page, eg. R10101
2. Click on the most specific pathway, in this case rn00760
3. Search the pathway diagram and find the reaction (highlighted in red). The arrows on the diagram indicate whether the reaction is reversible or irreversible, and which direction it goes.

I don't know how reliable these reversibilities are, nor the conditions in which they apply, so they should be treated with caution. However, as far as I know the reversibility is manually curated.

tl;dr painfully and not necessarily reliably.

[draeger]

The specification of KEGG's internally used representation format (KGML) says

The type attribute specifies the distinction of reversible and irreversible reactions, which are indicated by bi-directional and uni-directional arrows in the KEGG pathways. Note that the terms "reversible" and "irreversible" do not necessarily reflect biochemical properties of each reaction. They rather indicate the direction of the reaction drawn on the pathway map that is extracted from text books and literatures.

(see http://www.kegg.jp/kegg/xml/docs/)

[cdanielmachado]

@wbryant Thanks, I was not aware that the reversibility was indicated in the pathway maps. If it is indeed manually curated then I wonder how it was estimated and why it is not described in the reaction page itself.

Anyway, eQuilibrator gives a delta G estimate, you still need some assumptions to decide the reversibility from that. Could that be the reason for the inconsistency between your results and KEGG ?

[wbryant]

Thanks @draeger I hadn't come across the description of reversibility assignments. So given the breadth of sources for pathway maps the whole KEGG reversibility might not be internally consistent.

@cdanielmachado the delta G threshold could be the reason, I will have to look back over my code to check what assumptions I made.

[zakandrewking]

How useful would it be to have links from BiGG reaction & metabolite pages to Equilibrator?

[pstjohn]

That could definitely be useful! Probably wouldn't be too hard I'm assuming since the kegg links are already there.

pinging @jonstrutz11 who mentioned he was possibly working on a thermo calculator

[djinnome]

Also, the MetaCyc reaction & metabolite pages that are linked to from BiGG contain delta G calculations using the group contribution method.

[cdanielmachado]

The problem with BiGG to KEGG mappings is that they are not unique. I believe they come from MetaNetX (@zakandrewking can correct me if I'm wrong). You can often have one BiGG compound associated with multiple KEGG compounds, or no KEGG compound at all. So there is a lot of ambiguity (and loss of information) if one uses KEGG as an intermediate step to map between BiGG and eQuilibrator.

[djinnome]

Hi @cdanielmachado I can confirm that the mappings come from MetaNetX, which at least has a published methodology for performing the mappings. You are correct that KEGG Mappings are not unique or complete, but that is mainly due to the fact that BiGG metabolites do not have an InChI associated with them, so their true chemical structure is undefined. It must be inferred from the biochemical reactions in which it participates or from the IUPAC name.

[djinnome]

With that said, Adam Feist, Chris Henry, Markus Krummenaker and I spent quite a bit of time mapping the BiGG identifiers to EcoCyc and KEGG in order to compute deltaG0, and a published mol file is available for almost all the compounds in iAF1260 in supplementary information 5 that could be used to calculate the Gibbs free energy of formation for each metabolite using equilibrator or any other computational method.

[zakandrewking]

I'm tagging @bdu91 here. He might help out with linking to equilibrator.

[joanarcx]

Just want to add, this is of paramount importance to us here (Dusseldorf) too. Weren't Bin Du and Daniel Zielinski working in UCSD on an alternative tool to equilibrator for delta G estimation (they advertised it in the last CPH conference as both accurate and with high coverage:))? If that could go raw from BiGG, it would not only avoid problems with the mappings but give us hope in better thermodynamical estimations (I think this is the largest problem, no? equilibrator is still highly incomplete and unreliable for too many reactions).

[matthiaskoenig]

Using an annotated parameter to store Delta G values would allow to use the Delta G values in constraints, for instance to restrict the forward and backward fluxes of a reaction.

This could than work natively together with constraints in SBML fbc.

For instance via constraints like vf >= Delta_G * vb (does not make to have sense, just showing what one could do)

Via an annotation such things would be not possible to cleanly express. M

[Midnighter]

Some time has passed since the last comment but I wanted to update everyone here with the news that Elad and I reworked equilibrator to use all MetaNetX chemicals with a structure to build the database and secondly, that I'm working on functions to automatically map cobrapy metabolites and reactions to the equilibrator equivalent. This should enable plenty of other convenience functions.

Would love to hear what you guys think and in particular, your opinion on the preferred order of annotations to use for this.

• Merge request on equilibrator-cache that maps metabolites
• Merge request on equilibrator-api that translates reactions

[djinnome]

Much appreciated!

[cdanielmachado]

@Midnighter thanks for the update, this is great news. It would be extremely useful to be able to automatically get the delta G for any reaction in the BiGG database.

By the way, will this work with cross-compartment reactions?

[Midnighter]

By the way, will this work with cross-compartment reactions?

Probably best if @eladnoor responds himself but as far as I know, he is working on how to address transport reactions at the moment. So there is no special code to facilitate working with transports in equilibrator yet.

[eladnoor]

@Midnighter is right. We still don't have the implementation for it transport reactions. It is more complicated because one needs to define the aqueous parameters for both of the compartments. Since this data is not usually in the BiGG database or COBRA models, we need to get it from the user at the same time as the reactions themselves.

[rmtfleming]

Calculating thermochemical parameters, taking into account multiple compartments, can be achieved with the COBRA toolbox, using: https://opencobra.github.io/cobratoolbox/stable/tutorials/tutorialVonBertalanffy.html

To the best of my knowledge, under the hood of equilibrator and the COBRA Toolbox is the same component contribution method: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003098

The theory underlying Quantitative assignment of reaction directionality in a multicompartmental human metabolic reconstruction is here: https://www.ncbi.nlm.nih.gov/pubmed/22768925

On Wed, 11 Dec 2019 at 11:05, Elad Noor notifications@github.com wrote:

@Midnighter is right. We still don't have the implementation for it transport reactions. It is more complicated because one needs to define the aqueous parameters for both of the compartments. Since this data is not usually in the BiGG database or COBRA models, we need to get it from the user at the same time as the reactions themselves.

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

--

Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.

Assistant Professor, Division of Systems Biomedicine and Pharmacology, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University. https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming & H2020 Project Coordinator, Systems Medicine of Mitochondrial Parkinson’s Disease, http://sysmedpd.eu & Adjunct Lecturer, School of Medicine, National University of Ireland, Galway.

Peer-reviewed publications: https://goo.gl/FZPG23 Mobile: +353 873 413 072 Skype: ronan.fleming

Upcoming course from 25-29 November 2019, at Leiden University on "Constraint-based modelling: introduction and advanced topics". https://www.biosb.nl/archive-courses/constraint-based-modelling-introduction-and-advanced-topics-2019/

(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

[cdanielmachado]

I guess this is where the magic happens:

% Compartment specific pH
ph = [7.2; 7.4; 6.35; 5.5; 8; 7.2; 7.2; 7; 7.2];

% Compartment specific ionic strength in mol/L
is = 0.15*ones(length(compartments),1);

% Compartment specific electrical potential relative to cytosol in mV
chi = [0; 30; 0; 19; -155; 0; 0; -2.303*8.3144621e-3*T*(ph(compartments == 'x') - ph(compartments == 'c'))/(96485.3365e-6); 0];

Could this be easily generalized to the BiGG universal reaction database?

[eladnoor]

Calculating thermochemical parameters, taking into account multiple compartments, can be achieved with the COBRA toolbox

We were indeed planning in equilibrator to use the same framework as @rmtfleming suggests (i.e. vonBertalanffy, potentially ported to python). If the compartment information is added to COBRA models themselves, that would make things much easier. However, we are not sure how to determine the exact charge of the species that are transported (which is needed for the electrostatic term). How is this done in vonBertalanffy?

Could this be easily generalized to the BiGG universal reaction database?

I think that the problem is that the BiGG universal reaction database contains reactions without a specific organism/compartment context. I am not sure how we can use this database for ΔG'0 calculations without having a specific model (and all relevant parameters).

[eladnoor]

Calculating thermochemical parameters, taking into account multiple compartments, can be achieved with the COBRA toolbox...

I think @rmtfleming was referring to this part of the code: vonBertalanffy/estimateDrGt0.m. I do have a question about it though. It seems that the number of H atoms is taken directly from the formulas in the COBRA model. Can we be sure that they are set correctly? They are pH dependent, and one usually needs to use pKa data to determine them, but I guess the vonBertalanffy toolbox is the only place where that can happen. Are they set somewhere else in another script?

The theory underlying Quantitative assignment of reaction directionality in a multicompartmental human metabolic reconstruction...

And just one more question about the implementation. In equation 5 in this paper, Q refers to the "net number of charges transported from initial to final compartment", which I think happens in this line:

delta_chi = model.S'*spdiags(model.metCharges,0,nMet,nMet)*metCHI;

However, if the charge on both sides is not the same, this will not be equal to F*Q*ΔΨ, but to F*(Q2*Ψ2-Q1*Ψ1). How can we choose a single charge for the transport reaction in such cases? Or maybe that's not necessary?

[rmtfleming]

Hi Elad,

for transport reactions, it is assumed that the reaction given in the reconstruction is in terms of the species transported. Since the time of that paper, I have wondered how to systematically identify the species transported. Although one can estimate the pKa's and thereby predict the distribution of species in different protonated states, I know of no general biochemistry result that states that the predominant species is the one that transporters have evolved to transport. If the rate of transport of any protonated form of the same metabolite is slow with respect to the rate of equilibration of different protonated forms of the same metabolite, then transport could be selective for any one of a number of different species.

model.chi is a c x 1 array of compartment-specific electrical potential values in mV. model.metCharges is an m x 1 numerical array of metabolite charges. Both of these are expected to be provided with the reconstruction. i.e. we assume the reconstruction specifies the transported species. Moreover, we assume the reconstruction is charge-balanced.

I think that this line delta_chi = model.S'spdiags(model.metCharges,0,nMet,nMet)metCHI; is consistent with equation 29, derived in the supplementary material to the same paper.

In the supplementary material, Q is first defined as the charge of a species (model.metCharges in code). Later Q is defined as the net number of charges transported from initial to final compartment, but we did not state explicitly in that paragraph that we assumed that the reaction was charge-balanced, which I think we should have. However, in the main text, we state that "It is important to note that the stoichiometry of each transport reaction must accurately reflect the mechanism of intercompartmental transport, down to the details of the charged species transported." That means we expect the specification of the transport reaction to be charge-balanced. The contribution to the change in Gibbs energy comes from net movement of charge from one compartment to another, even though the reaction itself is charge balanced.

I hope that helps to clarify.

Regards,

Ronan

On Sat, 11 Jan 2020 at 11:40, Elad Noor notifications@github.com wrote:

Calculating thermochemical parameters, taking into account multiple compartments, can be achieved with the COBRA toolbox...

I think Ronan was referring to this part of the code: vonBertalanffy/estimateDrGt0.m. I do have a question about it though. It seems that the number of H atoms is taken directly from the formulas in the COBRA model. Can we be sure that they are set correctly? They are pH dependent, and one usually needs to use pKa data to determine them, but I guess the vonBertalanffy toolbox is the only place where that can happen. Are they set somewhere else in another script?

The theory underlying Quantitative assignment of reaction directionality in a multicompartmental human metabolic reconstruction...

And just one more question about the implementation. In equation 5 in this paper, Q refers to the "net number of charges transported from initial to final compartment", which I think happens in this line:

delta_chi = model.S'spdiags(model.metCharges,0,nMet,nMet)metCHI;

However, if the charge on both sides is not the same, this will not be equal to FQΔΨ, but to F(Q2Ψ2-Q1*Ψ1). How can we choose a single charge for the transport reaction in such cases? Or maybe that's not necessary?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

--

Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.

Assistant Professor, Division of Systems Biomedicine and Pharmacology, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University. https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming & H2020 Project Coordinator, Systems Medicine of Mitochondrial Parkinson’s Disease, http://sysmedpd.eu & Adjunct Lecturer, School of Medicine, National University of Ireland, Galway.

Peer-reviewed publications: https://goo.gl/FZPG23 Mobile: +353 873 413 072 Skype: ronan.fleming

(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

[eladnoor]

This makes things much clearer. One thing that's not completely clear to me yet, is how we can find out the transported species from the reconstruction itself? For example, say malate is transported from a compartment where its charge is -2 to one where it is -1. Since the reaction has to be charge balanced, the only way is to add a proton to the -2 side. That means that the total charge transported is -1. However, I don't see how I can change it to -2 (assuming I have specific information that is indeed the case). Perhaps I misunderstood how charge balancing is achieved?

[rmtfleming]

Hi Elad,

a malate example is here: https://www.vmh.life/#reaction/r2387 It is maybe not quite what you were looking for with your example.

Would a solution not be to associate with a proton in the new compartment?

met[x] + proton [y] <-> met[y]

Regards,

Ronan

On Wed, 15 Jan 2020 at 14:13, Elad Noor notifications@github.com wrote:

This makes things much clearer. One thing that's not completely clear to me yet, is how we can find out the transported species from the reconstruction itself? For example, say malate is transported from a compartment where its charge is -2 to one where it is -1. Since the reaction has to be charge balanced, the only way is to add a proton to the -2 side. That means that the total charge transported is -1. However, I don't see how I can change it to -2 (assuming I have specific information that is indeed the case). Perhaps I misunderstood how charge balancing is achieved?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/SBRG/bigg_models/issues/221?email_source=notifications&email_token=AAQMEOUYLABQQFIFAFTBUI3Q54K2JA5CNFSM4CYC5C22YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEJAOBVQ#issuecomment-574677206, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQMEOT7GKJBKC4SBMKCZRLQ54K2JANCNFSM4CYC5C2Q .

--

Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.

Assistant Professor, Division of Systems Biomedicine and Pharmacology, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University. https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming & H2020 Project Coordinator, Systems Medicine of Mitochondrial Parkinson’s Disease, http://sysmedpd.eu & Adjunct Lecturer, School of Medicine, National University of Ireland, Galway.

Peer-reviewed publications: https://goo.gl/FZPG23 Mobile: +353 873 413 072 Skype: ronan.fleming

(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

[eladnoor]

I am not sure if I got it right, but following your example, adding a proton will indeed add +1 the transported charge. What I was asking is if there is a way to make it more negative (we cannot really add a negative proton). By the way, another complication is that adding protons will also affect the reaction energetics if the ΔpH is not zero, wouldn't it?

[rmtfleming]

Hi Elad,

In this example, met[x], which is transported across a boundary, could have a charge of -2 and met[y] a charge of -1. met[x] + proton [y] <-> met[y]

I am not saying that this example above is always the way that a metabolite is transported if it has two different charges in two different compartments. An alternative would be: met[x] + proton [x] <-> met[y] where the transporter is selective for a metabolite with charge -1.

I think the question to ask is, what might be a systematic approach to the identification of the selectivity of a transporter for a particular charged state? Perhaps a molecular biophysics colleague of yours might be able to answer that. I do not have an answer.

Regards,

Ronan

On Wed, 15 Jan 2020 at 15:54, Elad Noor notifications@github.com wrote:

I am not sure if I got it right, but following your example, adding a proton will indeed add +1 the transported charge. What I was asking is if there is a way to make it more negative (we cannot really add a negative proton). By the way, another complication is that adding protons will also affect the reaction energetics if the ΔpH is not zero, wouldn't it?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

--

Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.

Assistant Professor, Division of Systems Biomedicine and Pharmacology, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University. https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming & H2020 Project Coordinator, Systems Medicine of Mitochondrial Parkinson’s Disease, http://sysmedpd.eu & Adjunct Lecturer, School of Medicine, National University of Ireland, Galway.

Peer-reviewed publications: https://goo.gl/FZPG23 Mobile: +353 873 413 072 Skype: ronan.fleming

(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

[eladnoor]

Sorry for my confusion. I didn't notice that you marked the substrate proton as the opposite compartment (y). That would solve the issue I raised, I think.

Regarding the broader issue of how to systematically assign the charges, I don't have any ideas either. I just wanted to first understand how one can encode it into the COBRA model.

Thanks for all the help!