In cdiener/dycone: dycone - analyze the dynamic landscape of metabolome data

Human blood cell example

In this example we will use dycone to analyze some perturbations in the Red blood cell model of Resendis. First we will use dycone to get its stochiometric matrix and extract the steady state from the orginial SBML file.

library(dycone)
data(eryth)
reacts = make_irreversible(eryth)
print(eryth)

spec = species(eryth)
concs = runif(length(spec))
names(concs) = spec
print(concs)

We extract the stochiometric matrix with the constant species simply by:

S = stoichiometry(eryth)

Extracting the k-cone

Getting the consistent k-cone (all $k_i>0$) takes two steps in dycone. First, calculation of the general flux cone and, second, calculation of the specific flux cone. Let's calculate and look at the flux cone first...

V = polytope_basis(S)
plot_basis(V)

Some interesting structurs can be observed from the flux kone alone. For instance, basis vectors with only two reactions indicate a constant relation of the form $k_i/k_j = const.$.

Perturbance experiment

Now, let's set up a small perturbation experiment in order to execute the differential k-cone analysis. We will use a faster version of the eigendynamics as reference parameters and perturb the steady state value of pyruvate.

m = ma_terms(S, concs)
K = kcone(V, m)
k = 1e3*eigendynamics(K)

# Now we generate several new steady states coming from a disturbed ethanol
# production

scales = 10^c(-1,-2)
d_concs = concs
for( i in 1:length(scales) ) {
    d_c = concs
    d_c["pyr"] = d_c["pyr"]*scales[i]
    d_concs = rbind(d_concs, d_c)
}

m_terms = apply(d_concs, 1, function(co) ma_terms(S, co))
d_basis = lapply(1:3, function(i) kcone(V, m_terms[,i]))

First we will have a look at the projection of the respective k-cones in a 2D space with colors from green to red indicating a diminishing ethanol production.

plot_red(d_basis)

This gives a visual representation but no detailed information which reactions are the most likely to be perturbed in the transition from the high to low pyruvate state. This can be achieved by the hyp function. We will add noise samples to each condition to simulate an experimental measurement.

noise = function(x) {
    new = pmax(0, x*rnorm(length(x),1,0.25))
    names(new) = names(x)
    return(new)
}

set.seed(42)

exp_concs = cbind(noise(d_concs[1,]), noise)
mats <- cbind(replicate(4, ma_terms(S, noise(d_concs[1,]))),
              replicate(4, ma_terms(S, noise(d_concs[3,]))))
samples <- factor(rep(c("normal", "disease"), each = 4))

p = minimal_perturbation(mats, samples, reacts)
print(p$k[p$k$adj.P.Val < 0.05,])

We see that the most prominent change is an increased reaction rate for two reactions, which result to be the two reactions consuming Pyruvate. However, there are also some other changes, however, their fold_changes are much smaller than 2, indicating only slight alterations due to the added noise. The mean log-fold change of the alterations in pyruvate production are close to the actual alteration (96-fold predicted vs 100-fold truth).

As we can see the change in Pyruvate is explained by either up-regulating the pyruvate consumption or (less likely) a lower production of upstream substrates which might have to act together. As we can see the reaction producing pyruvate is not in the list. This makes sense since it is a unidirectional reaction whose down-regulation would affect the steady state quite severely in the concentration of adp and pep, but not pyruvate, thus making it an unlikely target when observing a lower pyruvate concentration alone.

cdiener/dycone documentation built on May 13, 2019, 2:41 p.m.