8 Perturbations to metabolic networks
One of the most common uses of a metabolic model is to predict what happens when something is taken away: a gene knocked out, an enzyme inhibited by a drug. Which reactions can the cell not do without? Which pairs of reactions back each other up, so that either can be lost but not both? Constraint-based models answer such questions quickly, at the scale of whole genomes, and these predictions are used to find drug targets in pathogens and to design strains for metabolic engineering. This chapter builds on Chapter 7 and uses the same E. coli core model.
Chapter 10, Perturbations to Metabolic Networks: knock-outs (§10.1), gene versus reaction deletions and gene–protein–reaction rules (§10.1.1), synthetic lethals and algorithms to find them, including Fast-SL (§10.2), and applications such as drug target identification (§10.6). Chapter 9 covers methods that predict the flux distribution after a perturbation, such as MoMA (§9.2).
8.1 Single deletions
To delete a reaction from a constraint-based model, we fix its flux to zero by setting both its lower and upper bounds to zero (book, Eq. 10.2), and solve the FBA problem again. If the maximum growth rate falls below some small fraction of the wild type, 1% here, the deletion is called lethal, and the reaction essential.
Experiments delete genes, not reactions, and the two do not map one to one. Each reaction in the model carries a gene–protein–reaction (GPR) rule, a Boolean expression of the genes needed for it. An enzyme made of several subunits needs all of them (AND); isozymes can stand in for one another (OR). A gene deletion removes exactly those reactions whose rule becomes false. COBRApy and the COBRA Toolbox evaluate these rules for us.
We delete each of the 74 internal reactions of the model in turn (leaving out exchange reactions and the biomass reaction), and then each of its 137 genes, in the standard aerobic glucose medium.
Needs the COBRA Toolbox; run initCobraToolbox once per session before this script.
%% Single reaction and gene deletions in the E. coli core model
% A reaction is deleted by fixing its flux to zero (lower and upper bounds
% both 0; book Eq. 10.2); the maximum growth rate is then found by FBA. A
% deletion is called lethal if growth falls below 1% of the wild type.
% Genes are deleted through the gene-protein-reaction (GPR) rules, which
% the COBRA Toolbox evaluates for us: a reaction is lost only if the
% deletion makes its rule false, so isozymes can stand in for one another.
% Requires: COBRA Toolbox (run initCobraToolbox once per session first).
model = readCbModel(fullfile(sbp.data_dir(), 'e_coli_core.mat'));
sol = optimizeCbModel(model, 'max');
wt = sol.f;
% Delete internal reactions only: no exchange reactions, no biomass
rxns = sort(model.rxns(~startsWith(model.rxns, {'EX_', 'Biomass'})));
%% Reaction deletions
ratio = zeros(numel(rxns), 1);
for k = 1:numel(rxns)
m = changeRxnBounds(model, rxns{k}, 0, 'b');
sol = optimizeCbModel(m, 'max');
if sol.stat == 1
ratio(k) = sol.f / wt;
end % infeasible: no growth, ratio stays 0
end
essential_rxns = rxns(ratio < 0.01);
n_reduced = sum(ratio >= 0.01 & ratio < 0.99);
%% Gene deletions (the GPR rules decide which reactions are lost)
grRatio = singleGeneDeletion(model, 'FBA');
grRatio(isnan(grRatio)) = 0; % infeasible: no growth
[essential_genes, idx] = sort(model.genes(grRatio < 0.01));
names = model.geneNames(grRatio < 0.01);
names = names(idx);
fprintf('Wild-type growth %.4f /h; %d reactions tested\n', wt, numel(rxns));
fprintf('%d essential reactions: %s\n', numel(essential_rxns), ...
strjoin(essential_rxns, ', '));
fprintf('%d reactions whose loss reduces growth, but not below 1%%\n', ...
n_reduced);
fprintf('%d essential genes, of %d: %s\n', numel(essential_genes), ...
numel(model.genes), strjoin(strcat(essential_genes, ' (', names, ...
')'), ', '));
%% Plot: growth after each reaction deletion, sorted
[values, order] = sort(ratio);
colors = repmat([0.122 0.467 0.706], numel(values), 1); % no effect
colors(values < 0.99, :) = repmat([1.000 0.498 0.055], ...
sum(values < 0.99), 1); % reduced
colors(values < 0.01, :) = repmat([0.839 0.153 0.157], ...
sum(values < 0.01), 1); % lethal
fig = figure('Position', [100 100 1000 400]);
b = bar(values, 'FaceColor', 'flat');
b.CData = colors;
hold on
lethal = find(values < 0.01); % zero-height bars: mark them
plot(lethal, 0.03 * ones(size(lethal)), 'v', 'Color', [0.839 0.153 0.157], ...
'MarkerFaceColor', [0.839 0.153 0.157], 'MarkerSize', 5);
set(gca, 'XTick', 1:numel(values), 'XTickLabel', rxns(order), ...
'XTickLabelRotation', 90, 'FontSize', 7);
ylabel('Growth relative to wild type');
title('Single reaction deletions, E. coli core model, aerobic glucose');
box off
sbp.save_figure(fig, mfilename('fullpath'), 'fig-single-deletions');
results.wild_type_growth = wt;
results.growth_ratio = cell2struct(num2cell(ratio), rxns, 1);
results.essential_reactions = essential_rxns';
results.essential_genes = essential_genes';
sbp.save_results(mfilename('fullpath'), results);
Needs COBRApy.
"""Single reaction and gene deletions in the E. coli core model.
A reaction is deleted by fixing its flux to zero (lower and upper bounds
both 0; book Eq. 10.2); the maximum growth rate is then found by FBA. A
deletion is called lethal if growth falls below 1% of the wild type.
Genes are deleted through the gene-protein-reaction (GPR) rules, which
COBRApy evaluates for us: a reaction is lost only if the deletion makes
its rule false, so isozymes can stand in for one another.
"""
import numpy as np
import matplotlib.pyplot as plt
import cobra
from cobra.flux_analysis import single_gene_deletion
import sbp
model = cobra.io.read_sbml_model(str(sbp.DATA / "e_coli_core.xml"))
wt = model.slim_optimize()
cutoff = 0.01 * wt
# Delete internal reactions only: no exchange reactions, no biomass
rxns = sorted(r.id for r in model.reactions
if not r.id.startswith(("EX_", "Biomass")))
# --- Reaction deletions
ratio = {}
for rid in rxns:
with model: # undone after the block
model.reactions.get_by_id(rid).bounds = (0, 0)
ratio[rid] = model.slim_optimize(error_value=0.0) / wt
essential_rxns = [r for r in rxns if ratio[r] < 0.01]
reduced = [r for r in rxns if 0.01 <= ratio[r] < 0.99]
# --- Gene deletions (the GPR rules decide which reactions are lost)
gd = single_gene_deletion(model)
gene_ratio = {next(iter(ids)): (g if np.isfinite(g) else 0.0) / wt
for ids, g in zip(gd["ids"], gd["growth"])}
essential_genes = sorted(g for g, r in gene_ratio.items() if r < 0.01)
names = [model.genes.get_by_id(g).name for g in essential_genes]
print(f"Wild-type growth {wt:.4f} /h; {len(rxns)} reactions tested")
print(f"{len(essential_rxns)} essential reactions: "
f"{', '.join(essential_rxns)}")
print(f"{len(reduced)} reactions whose loss reduces growth, but not below 1%")
print(f"{len(essential_genes)} essential genes, of {len(model.genes)}: "
f"{', '.join(f'{g} ({n})' for g, n in zip(essential_genes, names))}")
# --- Plot: growth after each reaction deletion, sorted
order = sorted(rxns, key=lambda r: ratio[r])
values = np.array([ratio[r] for r in order])
colors = np.where(values < 0.01, "#d62728",
np.where(values < 0.99, "#ff7f0e", "#1f77b4"))
sbp.style()
fig, ax = plt.subplots(figsize=(10, 4))
ax.bar(range(len(order)), values, color=colors)
lethal = np.flatnonzero(values < 0.01) # zero-height bars: mark them
ax.plot(lethal, np.full(lethal.size, 0.03), "v", color="#d62728", ms=5)
ax.set_xticks(range(len(order)), order, rotation=90, fontsize=7)
ax.set_ylabel("Growth relative to wild type")
ax.set_title("Single reaction deletions, E. coli core model, aerobic glucose")
for col, lab in [("#d62728", "lethal (< 1%)"), ("#ff7f0e", "reduced"),
("#1f77b4", "no effect")]:
ax.bar(0, 0, color=col, label=lab)
ax.legend(frameon=False, loc="upper left")
sbp.save_figure(fig, __file__, "fig-single-deletions")
sbp.save_results(__file__, {
"wild_type_growth": wt,
"growth_ratio": ratio,
"essential_reactions": essential_rxns,
"essential_genes": essential_genes,
})
sbp.show()
Thirteen reactions are essential (red), among them the glucose uptake system (GLCpts), the lower steps of glycolysis (GAPD, PGK, PGM, ENO), the first steps of the TCA cycle (CS, ACONTa/b, ICDHyr), and the uptake of ammonium and phosphate. Losing another 23 reactions slows growth without stopping it (orange): the network reroutes flux, at a cost. Deleting the ATP maintenance reaction, ATPM, even increases growth slightly, since the ATP it consumes is freed for biosynthesis.
Only seven genes are essential, although thirteen reactions are. The difference comes from the GPR rules. Aconitase (ACONTa and ACONTb), for example, is encoded by two genes, acnA and acnB, either of which suffices, so neither gene is essential although the reaction is. In the other direction, ptsH and ptsI are each essential because the glucose transport system needs both. Predictions about genes, which are what experiments test, can therefore differ substantially from predictions about reactions.
8.2 Synthetic lethals
Two reactions form a synthetic lethal pair if deleting either one alone leaves the cell viable, but deleting both is lethal (book, §10.2). Such pairs reveal how a network backs itself up, and they are attractive targets for combination drugs.
The obvious way to find them is to test every pair of non-essential reactions. For the core model that is \(\binom{61}{2} = 1830\) double deletions, which is easy; for a genome-scale model with 2000 reactions it is about two million, and for triplets more than a billion. The Fast-SL algorithm, developed in our group, cuts down the search with a simple argument (Pratapa et al. 2015). If some flux distribution that supports growth uses neither reaction of a pair, deleting both leaves that distribution intact, so the pair cannot be lethal. Hence:
- at least one reaction of a lethal pair, call it \(i\), must carry flux in a wild-type solution; and
- the other reaction, \(j\), must carry flux in a solution that remains once \(i\) is deleted.
The sparser these solutions, the fewer pairs need to be tested, so Fast-SL uses solutions that achieve the maximum growth with the minimum total flux, \(\sum_j |v_j|\) (parsimonious FBA; book, §8.8.5). The code below finds all synthetic lethal pairs both ways, and checks that the answers agree.
Needs the COBRA Toolbox.
%% Synthetic lethal reaction pairs in the E. coli core model
% A pair of reactions is synthetic lethal if deleting both stops growth
% (below 1% of wild type) while deleting either one alone does not. We
% find all such pairs twice:
% 1. exhaustively, testing every pair of non-essential reactions;
% 2. as in Fast-SL (Pratapa et al., 2015), testing only pairs that can
% possibly be lethal. If a growing flux distribution uses neither
% reaction of a pair, deleting both leaves it intact, so the pair is
% not lethal. Hence at least one reaction of a lethal pair carries flux
% in the wild-type solution, and the other carries flux once the first
% is deleted. Using sparse solutions (minimal total flux) keeps the
% number of pairs to test small.
% Requires: COBRA Toolbox (run initCobraToolbox once per session first).
model = readCbModel(fullfile(sbp.data_dir(), 'e_coli_core.mat'));
sol = optimizeCbModel(model, 'max');
cutoff = 0.01 * sol.f;
rxns = sort(model.rxns(~startsWith(model.rxns, {'EX_', 'Biomass'})));
is_essential = false(size(rxns));
for k = 1:numel(rxns)
is_essential(k) = growth_without(model, rxns(k)) < cutoff;
end
candidates = rxns(~is_essential);
%% 1. Exhaustive: every pair of non-essential reactions
pairs = nchoosek(1:numel(candidates), 2);
lethal_exhaustive = {};
for p = 1:size(pairs, 1)
pair = candidates(pairs(p, :))';
if growth_without(model, pair) < cutoff
lethal_exhaustive(end + 1, :) = pair; %#ok<SAGROW>
end
end
%% 2. Fast-SL: only pairs (i, j) with i active in the wild type and j
% active once i is deleted
tested = {};
lethal_fastsl = {};
for i = active_without(model, {}, candidates)'
for j = active_without(model, i, candidates)'
pair = sort([i, j]);
key = strjoin(pair, ' ');
if ~strcmp(i{1}, j{1}) && ~ismember(key, tested)
tested{end + 1} = key; %#ok<SAGROW>
if growth_without(model, pair) < cutoff
lethal_fastsl(end + 1, :) = pair; %#ok<SAGROW>
end
end
end
end
lethal_fastsl = sortrows(lethal_fastsl);
same = isequal(lethal_fastsl, sortrows(lethal_exhaustive));
fprintf('%d non-essential reactions\n', numel(candidates));
fprintf('Exhaustive: %d pairs tested, %d synthetic lethal\n', ...
size(pairs, 1), size(lethal_exhaustive, 1));
fprintf('Fast-SL: %d pairs tested, %d synthetic lethal\n', ...
numel(tested), size(lethal_fastsl, 1));
fprintf('Same pairs found: %d\n', same);
partners = cellfun(@(r) sum(strcmp(lethal_fastsl(:), r)), candidates);
[~, order] = sortrows([-partners, (1:numel(candidates))']);
top = order(1:5);
fprintf('Reactions in the most pairs: %s\n', strjoin(strcat( ...
candidates(top), ' (', arrayfun(@num2str, partners(top), ...
'UniformOutput', false), ')'), ', '));
%% Plot: which reactions form lethal pairs
involved = unique(lethal_fastsl(:));
M = zeros(numel(involved));
for p = 1:size(lethal_fastsl, 1)
a = strcmp(involved, lethal_fastsl{p, 1});
b = strcmp(involved, lethal_fastsl{p, 2});
M(a, b) = 1; M(b, a) = 1;
end
fig = figure('Position', [100 100 1100 500]);
subplot(1, 3, [1 2])
imagesc(M, [0 1]);
colormap([0.97 0.98 1; 0.12 0.47 0.71]);
axis square
set(gca, 'XTick', 1:numel(involved), 'XTickLabel', involved, ...
'XTickLabelRotation', 90, 'YTick', 1:numel(involved), ...
'YTickLabel', involved, 'FontSize', 7, 'Box', 'off');
title(sprintf('%d synthetic lethal pairs', size(lethal_fastsl, 1)));
subplot(1, 3, 3)
b = bar([size(pairs, 1), numel(tested)], 'FaceColor', 'flat');
b.CData = [0.6 0.6 0.6; 0.122 0.467 0.706];
set(gca, 'XTickLabel', {'exhaustive', 'Fast-SL'});
ylabel('Double deletions simulated');
title('Cost of the search');
box off
sbp.save_figure(fig, mfilename('fullpath'), 'fig-synthetic-lethals');
results.n_candidates = numel(candidates);
results.n_pairs_exhaustive = size(pairs, 1);
results.lethal_pairs = num2cell(lethal_fastsl, 2)';
results.same_as_exhaustive = same;
results.n_pairs_fastsl = numel(tested);
sbp.save_results(mfilename('fullpath'), results);
%% Local functions
function g = growth_without(model, deleted)
model = changeRxnBounds(model, deleted, 0, 'b');
sol = optimizeCbModel(model, 'max');
g = 0;
if sol.stat == 1
g = sol.f;
end
end
function active = active_without(model, deleted, candidates)
% Candidate reactions carrying flux in a minimal-total-flux solution
if ~isempty(deleted)
model = changeRxnBounds(model, deleted, 0, 'b');
end
sol = optimizeCbModel(model, 'max', 'one');
v = sol.x(findRxnIDs(model, candidates));
active = candidates(abs(v) > 1e-9);
end
Needs COBRApy.
"""Synthetic lethal reaction pairs in the E. coli core model.
A pair of reactions is synthetic lethal if deleting both stops growth
(below 1% of wild type) while deleting either one alone does not. We
find all such pairs twice:
1. exhaustively, testing every pair of non-essential reactions;
2. as in Fast-SL (Pratapa et al., 2015), testing only pairs that can
possibly be lethal. If a growing flux distribution uses neither
reaction of a pair, deleting both leaves it intact, so the pair is
not lethal. Hence at least one reaction of a lethal pair carries flux
in the wild-type solution, and the other carries flux once the first
is deleted. Using sparse solutions (minimal total flux, pFBA) keeps
the number of pairs to test small.
"""
import itertools
import numpy as np
import matplotlib.pyplot as plt
import cobra
from cobra.flux_analysis import pfba
import sbp
model = cobra.io.read_sbml_model(str(sbp.DATA / "e_coli_core.xml"))
wt = model.slim_optimize()
cutoff = 0.01 * wt
rxns = sorted(r.id for r in model.reactions
if not r.id.startswith(("EX_", "Biomass")))
def growth_without(deleted):
with model:
for rid in deleted:
model.reactions.get_by_id(rid).bounds = (0, 0)
return model.slim_optimize(error_value=0.0)
def active_without(deleted, candidates):
"""Candidate reactions carrying flux in a minimal-total-flux solution."""
with model:
for rid in deleted:
model.reactions.get_by_id(rid).bounds = (0, 0)
v = pfba(model).fluxes
return [r for r in candidates if abs(v[r]) > 1e-9]
essential = [r for r in rxns if growth_without([r]) < cutoff]
candidates = [r for r in rxns if r not in essential]
# --- 1. Exhaustive: every pair of non-essential reactions
pairs = list(itertools.combinations(candidates, 2))
lethal_exhaustive = [p for p in pairs if growth_without(p) < cutoff]
# --- 2. Fast-SL: only pairs (i, j) with i active in the wild type and j
# active once i is deleted
tested = set()
lethal_fastsl = set()
for i in active_without([], candidates):
for j in active_without([i], candidates):
pair = tuple(sorted((i, j)))
if j != i and pair not in tested:
tested.add(pair)
if growth_without(pair) < cutoff:
lethal_fastsl.add(pair)
lethal_fastsl = sorted(lethal_fastsl)
print(f"{len(candidates)} non-essential reactions")
print(f"Exhaustive: {len(pairs)} pairs tested, "
f"{len(lethal_exhaustive)} synthetic lethal")
print(f"Fast-SL: {len(tested)} pairs tested, "
f"{len(lethal_fastsl)} synthetic lethal")
print(f"Same pairs found: {lethal_fastsl == sorted(lethal_exhaustive)}")
partners = {r: sum(r in p for p in lethal_fastsl) for r in candidates}
top = sorted(partners, key=lambda r: (-partners[r], r))[:5]
print("Reactions in the most pairs: "
+ ", ".join(f"{r} ({partners[r]})" for r in top))
# --- Plot: which reactions form lethal pairs
involved = sorted({r for p in lethal_fastsl for r in p})
idx = {r: k for k, r in enumerate(involved)}
M = np.zeros((len(involved), len(involved)))
for a, b in lethal_fastsl:
M[idx[a], idx[b]] = M[idx[b], idx[a]] = 1
sbp.style()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 5),
gridspec_kw={"width_ratios": [2.2, 1]})
ax1.imshow(M, cmap="Blues", vmin=0, vmax=1.3)
ax1.set_xticks(range(len(involved)), involved, rotation=90, fontsize=7)
ax1.set_yticks(range(len(involved)), involved, fontsize=7)
ax1.set_title(f"{len(lethal_fastsl)} synthetic lethal pairs")
for spine in ax1.spines.values():
spine.set_visible(False)
ax2.bar(["exhaustive", "Fast-SL"], [len(pairs), len(tested)],
color=["#999999", "#1f77b4"])
ax2.set_ylabel("Double deletions simulated")
ax2.set_title("Cost of the search")
sbp.save_figure(fig, __file__, "fig-synthetic-lethals")
sbp.save_results(__file__, {
"n_candidates": len(candidates),
"n_pairs_exhaustive": len(pairs),
"lethal_pairs": [list(p) for p in lethal_fastsl],
"same_as_exhaustive": lethal_fastsl == sorted(lethal_exhaustive),
"n_pairs_fastsl": len(tested),
})
sbp.show()
Both searches find the same 88 synthetic lethal pairs, but Fast-SL simulates only about a quarter as many double deletions. The saving grows rapidly with the size of the network and with the order of the deletions, which is what makes searches for lethal triplets and quadruplets in genome-scale models feasible.
Many of the pairs make biological sense. The reactions of upper glycolysis (PGI, PFK, FBA, TPI) are each paired with reactions of the pentose phosphate pathway (G6PDH2r, PGL, GND, RPE, TALA, TKT1, TKT2): the two pathways are alternative routes from glucose 6-phosphate into lower glycolysis, and the cell can use either, but not neither. Glutamate can be made by glutamate dehydrogenase (GLUDy) or by glutamate synthase (GLUSy), and anaplerotic PEP carboxylase (PPC) is backed up by the glyoxylate shunt (ICL, MALS). Other pairs, such as those involving water transport (H2Ot), arise from the way the model is written, and are a reminder that every prediction must be checked against biology (book, §10.7).
The exact number of double deletions that Fast-SL needs depends on which of several equally sparse solutions the solver returns, so it can differ slightly between MATLAB and Python; the set of lethal pairs cannot, and the automated checks require it to be identical.
8.3 Exercises
- Repeat the single deletions for anaerobic growth (oxygen uptake set to zero). Which reactions become essential, and why?
- Use the gene-deletion functions (
single_gene_deletionin COBRApy,singleGeneDeletionin the COBRA Toolbox) to find essential genes, and the GPR rules to explain every case where a reaction is essential but none of its genes is. - Change the lethality cut-off from 1% to 5% and to 50% of wild-type growth. How do the lists of essential reactions and synthetic lethal pairs change?
- Extend the Fast-SL code to find synthetic lethal triplets: three reactions whose joint deletion is lethal, when no pair among them is. How many triple deletions does it test, compared with the \(\binom{n}{3}\) of an exhaustive search? (For genome-scale models, use the full implementations at https://github.com/RamanLab/FastSL and https://github.com/RamanLab/FastSL-py.)
- A minimal reactome is a set of reactions that supports growth, from which no reaction can be removed without killing the cell. Starting from all reactions, remove reactions one at a time, in random order, keeping each removal only if the cell still grows. How large is the result, and does it change with the order? (This is the second problem of the BT5240 constraint-based assignment.)
8.4 Going further
Deletions are the simplest perturbations. Over-expression, enzyme inhibition and changes of medium can be modelled in the same framework, and methods such as MoMA and ROOM predict how fluxes redistribute after a perturbation, instead of assuming that the mutant immediately grows optimally (book, Chapters 9 and 10).