Pythooooon

I’ve written previously about the Daaaaave algorithm for using ‘omics data to predict fluxes through a metabolic network. We originally implemented the algorithm in MatLab, and released the code as supplementary material to the paper [1]. But the journal (BMC Systems Biology) messed up: changing the file names, and breaking the code.

We’ve tried to get them to fix this. Unsuccessfully.

So instead, I’ve put the code on sourceforge github [2]. And since MatLab sits uncomfortably with my cardigan-wearing, open source approach to science, I’ve reimplemented the algorithm in python there, too.

References

1. Lee D, Smallbone K, Dunn WB, Murabito E, Winder CL, Kell DB, Mendes P, Swainston N (2012) “Improving metabolic flux predictions using absolute gene expression data” BMC Systems Biology 6:73.
   doi:10.1186/1752-0509-6-73
2. http://github.com/u003f/daaaaave

A python

Advertisement

β-oxing

Dietary fibres are fermented by bacteria in the colon to form short-chain fatty acids: acetate (C₂), propionate (C₃) and butyrate (C₄). Fatty acids are taken up by colonocytes, and processed in the mitochondria using β-oxidation. This metabolism of — otherwise indigestible — carbohydrates via gut flora provides a major source of energy in humans.

Those of you interested in SCFA metabolism would do well to follow the work of Barbara Bakker [1-3]. In a recent paper with Karen van Eunen and friends, Babs created a mathematical model of β-ox in the rat liver. The model considers competition between acyl-CoAs with differing chain lengths for the same enzymes, which can lead to pathway saturation. I’ve created SBML versions of their model: have a play.

• “steady-state model”: BIOMD0000000505
• “time-course model”: BIOMD0000000506

overview of β-ox [1]

References

1. van Eunen K, Simons SM, Gerding A, Bleeker A, den Besten G, Touw CM, Houten SM, Groen BK, Krab K, Reijngoud DJ, Bakker BM (2013) “Biochemical competition makes fatty-acid β-oxidation vulnerable to substrate overload” PLoS Computational Biology 9:e1003186.
   doi:10.1371/journal.pcbi.1003186
2. den Besten G, van Eunen K, Groen AK, Venema K, Reijngoud DJ, Bakker BM (2013) “The role of short-chain fatty acids in the interplay between diet, gut microbiota, and host energy metabolism” Journal of Lipid Research 54:2325-2340.
   doi:10.1194/jlr.R036012
3. den Besten G, Lange K, Havinga R, van Dijk TH, Gerding A, van Eunen K, Müller M, Groen AK, Hooiveld GJ, Bakker BM, Reijngoud DJ (2013) “Gut-derived short-chain fatty acids are vividly assimilated into host carbohydrates and lipids” American Journal of Physiology – Gastrointestinal and Liver Physiology 305:G900-910.
   doi:10.1152/ajpgi.00265.2013

A Baltic Circuit

Last week I took the train and ferry to Sønderborg, in Denmark. Whilst there, U+003F junior went to his first gig; it was Efterklang’s last.

I like trains. I like the Baltics. I like ferries and Efterklang too, but let’s not complicate matters. You can [almost] circumnavigate the Baltic Sea by train. Here’s how:

• Day 1: starting at my little brother’s place in Göteborg, take the sleeper at 18:25, arriving into Luleå in the north of Sweden at 11:30 the next day. Cost: 873kr.
• Day 2: take the bus [boo!] across the border into Finland, leaving at 12:45 and arriving to Tornio at 15:15. Cost: free! Then a sleeper to Helsinki, leaving at 20:31 and arriving at 10:00. Cost: €105.
• Day 3: travel to Moscow on the Tolstoi, leaving at 17:23, arriving into Moscow Leningradski at 08:24. This is €158, but you’ll also need to get a Russian visa well in advance (£116). Spend the night in Moscow.
• Day 5: take the Trans-European Express to Berlin (18000₽), leaving Moscow Byelorruski at 07:44, arriving 06:53. This one runs through Belarus, so you’ll need another visa (£50).
• Day 6: The Berlin Night Express is a unique train that travels aboard a ferry to cross the Baltic. 799kr. Leaves Berlin at 22:30, arrives into Malmö at 12:55. It’s a short hop back to Göteborg, leaving at 13:08, arriving 16:20 on day 7. 343kr.

Total time: just under six days. Total cost: just over £900, based on two sharing etc etc, plus a Moscow hotel. To get back where you started. Brilliant.

The Baltic Sea

From genes to fluxes #2

@u003f Good point. Note though that in (Becker, 2008) the minimum and the maximum are used.

—
Nikos Ignatiadis (@nikosIgnatiadis) February 13, 2014

Last time I talked about the issue with the method used by Daaaaave [1] to map from enzyme-level data to reaction-level data. Given enzyme (or gene) levels A = 4, B = 3 and C = 2 units, we find:

reaction	GPR	Daaaaave	Gimme
4	(A and B) or (A and C)	min(4,3) + min(4,2) = 5	max(min(4,3), min(4,2)) = 3
5	A and (B or C)	min(4,3 + 2) = 4	min(4, max(3,2)) = 3

.

The problem with applying the min/plus rule to GPRs is that reactions 4 and 5 are the same (albeit differently bracketed), but Daaaaave assigns them different values. As Nikos pointed out, the min/max rule used by Gimme [2] doesn’t make this mistake. However, I think we really should be adding the activities of alternative catalysts; indeed some networks — such as “Yeast 1” [3] — use separate reactions in place of “or” statements. Any mapping must be robust to equivalent representations.

Let’s step back a bit. Reaction 4 is catalysed by alternative complexes, A:B and A:C.

r₄ → (A and B) or (A and C)

There is less of A (4) than the total amount of B (3) and C (2), so there must be some B or C “wasted” when forming the two complexes. There are an infinite number of arrangements here — we could have A:B/A:C = 3/1, 2/2, 2½/1½, … — but their maximum total activity is 4 units. This value of 4 is overestimated by Daaaaave, but underestimated by Gimme.

We can frame our verbal reasoning above mathematically. Each GPR mention of an enzyme across the network is really a separate entity

r₄ → (A₁ and B₁) or (A₂ and C₁)

that together make up the total enzyme level

A₁ + A₂ + … = A = 4.

We can substitute “and” relationships by introducing new variables X_i ≥ 0 that represent complexes

r₄ → X₁ or X₂

whose activities can be no more than any of their parts

X₁ ≤ A₁, X₁ ≤ B₁.

We can also substitute “or” relationships by introducing new variables Y_i that represent alternative catalysts

r₄ = Y₁

whose activities are the sum of their parts

Y₁ = X₁ + X₂.

Finally, we want there to be as little wastage as possible, and one way to achieve this is through maximising the total activity

maximise: r₁ + r₂ + ….

This optimisation is an LP problem and can be easily solved for networks of any size. Indeed, running an FBA over the network would be of the same computational complexity. Most importantly, this mapping makes the most of the available data.

References

1. Lee D, Smallbone K, Dunn WB, Murabito E, Winder CL, Kell DB, Mendes P, Swainston N (2012) “Improving metabolic flux predictions using absolute gene expression data” BMC Systems Biology 6:73.
   doi:10.1186/1752-0509-6-73
2. Becker SA, Palsson BØ (2008) “Context-specific metabolic networks are consistent with experiments” PLoS Comp Biol 4:e1000082.
   doi:10.1371/journal.pcbi.1000082
3. Herrgård MJ, Swainston N, Dobson P, Dunn WB, Arga KY, Arvas M, Blüthgen N, Borger S, Costenoble R, Heinemann M, Hucka M, Le Novère N, Li P, Liebermeister W, Mo ML, Oliveira AP, Petranovic D, Pettifer S, Simeonidis E, Smallbone K, Spasić I, Weichart D, Brent R, Broomhead DS, Westerhoff HV, Kirdar B, Penttilä M, Klipp E, Palsson BØ, Sauer U, Oliver SG, Mendes P, Nielsen J, Kell DB (2008) “A consensus yeast metabolic network reconstruction obtained from a community approach to systems biology” Nat Biotechnol 26:1155-1160.
   doi:10.1038/nbt1492

From genes to fluxes

In 2012, we set out a pipeline [1] for using predicting flux through a metabolic network, using genetic or proteomic data.

Two inputs are required:

• a genome-scale metabolic network containing gene-protein-reaction (GPR) associations, such as human or yeast
• quantitative estimates of their enzyme levels, via proteomics or absolute transcriptomic data

To predict metabolic flux, two steps are taken:

• enzyme-level data are combined with GPR associations, to create reaction-level data
• reaction-level data are constrained by mass-balancing, to create system-level flux data

A number of methods have been proposed to perform this second mapping, maximising the correlation between reaction data and system flux, including Daaaaave [1, developed in Manchester], Gimme [2] and iMat [3]. I won’t discuss their relative merits here, but will instead focus on the first step.

The mapping isn’t quite as easy as it seems, due to the many-to-many relationship between enzymes and metabolic reactions. Suppose we have determined enzyme levels of A = 4, B = 3 and C = 2. How should we interpret their various combinations as GPRs?

reaction	GPR	level
1	A	4
2	A or B	4 + 3 = 7
3	A and B	min(4,3) = 3
4	(A and B) or (A and C)	min(4,3) + min(4,2) = 5
5	A and (B or C)	min(4,3 + 2) = 4

.

Given the simple mapping in reaction 1, we may use the enzyme level directly. The “or” relationship in reaction 2 allows for alternative catalysts, so its total level is given by the sum of its components. The “and” relationship in reaction 3 means it is catalysed by a complex, whose maximum possible concentration is given by the minimum level of its components. These two min/plus rules may be combined for more complex GPRs such as in reactions 4 and 5.

Here’s the rub: the Boolean logic used in 4 and 5 are the same (just bracketed differently), yet the levels output by the mapping are different. The existing approach is badly-defined and we need to go back to the drawing-board. One option — proposed by Brandon Barker — is to insist that GPRs are always written in a consistent manner. For example, we could use the disjunctive normal form (DNF), meaning that GPRs are expressed as an “or of ands”: a list of alternative complexes, as per reaction 4. I think that having consistency like this is an excellent idea. But I also think that our mapping should be robust to alternative bracketing; I’ll show you how, next time.

References

1. Lee D, Smallbone K, Dunn WB, Murabito E, Winder CL, Kell DB, Mendes P, Swainston N (2012) “Improving metabolic flux predictions using absolute gene expression data” BMC Systems Biology 6:73.
   doi:10.1186/1752-0509-6-73
2. Becker SA, Palsson BØ (2008) “Context-specific metabolic networks are consistent with experiments” PLoS Comp Biol 4:e1000082.
   doi:10.1371/journal.pcbi.1000082
3. Shlomi T, Cabili MN, Herrgård MJ, Palsson BØ, Ruppin E (2008) “Network-based prediction of human tissue-specific metabolism” Nature Biotechnology 26:1003–1010.
   doi:10.1038/nbt.1487

Striking a balance with Recon 2.1

In 2013, after years of work, a comprehensive reconstruction of human metabolism, called “Recon 2“, was released to much fanfare [1].

But, as with any many-cooked project, some errors slipped through. The biggest issue I encountered involves unbalanced reactions and generic metabolites. Here are three reactions from the model:

1. generic fatty acid R-group transfer from CoA (C₂₁) to DHAP (C₃)

This reaction provides a template to be followed by all fatty acid side groups. However, to function as a model, there need to be connections between these generic molecules and their specific counterparts:

2. specific to generic reaction, converting C₁₄–CoA to R–CoA
3. generic-to-specific reaction, converting R–CoA to C₁₆–CoA

Taken together, reactions 2 and 3 are elementally-unbalanced, providing a mechanism for generating two carbons for free, and meaning any analyses with Recon 2 examining carbon flux  will be unreliable.

We can, however, modify the network to ensure that all reactions are elementally-balanced. Generic molecules must be replaced by their specific counterparts, creating a new reaction for each possible R-group, and removing any specific-to-generic interconversions:

Given n possible R-groups, there are n³ possible triglycerides. Replacing all generic metabolites as above thus leads to a much larger network, “Recon 2.1”.

	Recon 2	Recon 2.1
reactions	7440	71159
variables	5463	13940

.

The number of unique metabolites still falls well short of the forty thousand defined in HMDB, but is a step in the right direction.

The new model is available in SBML format from BioModels.net (MODEL1311110001). A more detailed description (and further tinkerings) can be found on arXiv [2].

References

1. Thiele I, Swainston N, Fleming RM, Hoppe A, Sahoo S, Aurich MK, Haraldsdottir H, Mo ML, Rolfsson O, Stobbe MD, Thorleifsson SG, Agren R, Bölling C, Bordel S, Chavali AK, Dobson P, Dunn WB, Endler L, Hala D, Hucka M, Hull D, Jameson D, Jamshidi N, Jonsson JJ, Juty N, Keating S, Nookaew I, Le Novère N, Malys N, Mazein A, Papin JA, Price ND, Selkov E Sr, Sigurdsson MI, Simeonidis E, Sonnenschein N, Smallbone K, Sorokin A, van Beek JH, Weichart D, Goryanin I, Nielsen J, Westerhoff HV, Kell DB, Mendes P, Every Man, His Dog, Palsson BØ (2013) “A community-driven global reconstruction of human metabolism” Nature Biotechnology 31:419-425.
   doi:10.1038/nbt.2488
2. Smallbone K (2013) “Striking a balance with Recon 2.1” arXiv:1311.5696

Surreal goings-on at MMU

Surreal numbers are a construction developed by British mathematician John Conway (who also brought us the Game of Life). They are defined recursively: given two surreal numbers x_L less than x_R, we can define a new number X = (x_L|x_R) that bisects the two. Starting with the empty set and using this rule, we can generate all the real numbers, and more.

But, of course, those of you from Manchester knew all this already, as their definition may (surreally) be found on Chester Road, at the entrance to the MMU car park.

This final panel is my favourite. It seems to say that starting from nothing (∅) we can produce everything (∞),  thus proving the existence a higher power (JHWH, the Tetragrammaton). In fact it is being deliberately ambiguous; “JHWH” is also how American computer scientist Donald Knuth referred to Conway, writing in 1979:

In the beginning, everything was void, and J.H.W.H.Conway began to create numbers. Conway said, “Let there be two rules which bring forth all numbers large and small. This shall be the first rule: Every number corresponds to two sets of previously created numbers, such that no member of the left set is greater than or equal to any member of the right set. And the second rule shall be this: One number is less than or equal to another number if and only if no member of the first number’s left set is greater than or equal to the second number, and no member of the second number’s right set is less than or equal to the first number.” And Conway examined these two rules he had made, and behold! they were very good.

The mystery remains as to how this mathematical memorial found its way to such an unexpected location. Answers on a postcard.

Thanks to my collaborator and lunch-mate Jon Borresen for pointing out the gates.

In the beginning there was #postaday

The days of heading to your University library to pore over the latest issue of Acta Arithmetica are a thing of the past. It’s a shame, but Generation Z want their research delivered electronically, and they want it delivered immediately. They don’t have the patience for polishing, for peer review, or for paper. For the academic, blogging has a role to play in this new world order, finding a niche in the immediacy:detail spectrum somewhere between twitter (one sentence critiques) and arXiv (unreviewed reports).

I first entered the blogoshpere in 2010. Inspired by gepasi‘s photo-a-day project 365, I embarked on a blogging #postaday challenge. Of course I failed. To write a decent blog post every day, you’d have to be highly imaginative, highly driven, and highly unemployed. I remain none of the above. But I’m back, armed with four years’ material and an 86% less demanding task. Wish me luck!

#postaweek