From b02f80f7c9ac085d87b9eb4d251b187a963ee368 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 22 Jun 2026 02:21:36 +0200 Subject: [PATCH] Add assignCompartments: functionality-constrained compartment-assignment MILP MATLAB port of raven-toolbox/assignCompartments' core: a single MILP (built as a COBRA prob and solved via optimizeProb, like gapFillMILP) that places the requested reactions into the compartments named in a GSS while requiring the objective to stay producible. Pathway coherence is emergent -- a reaction is placed against its own top score when the network needs it there. Mono-localization; passive transports (default<->c) added where needed; merges the model to a single compartment first (like predictLocalization), then re-distributes. optimizeProb already sets intTol=1e-9, so the Big-M flux gating does not leak. Returns the compartmentalised model, the per-reaction placement, and the added transports. Validated in MATLAB (tAssignCompartments, 3 tests): functionality overrides the localization score (a confined metabolite keeps a reaction in cytosol against its higher score), the score is respected when transports keep the model functional, and an unreachable growth floor is reported. Gap-fill coupling (universal/z candidates) is deferred to a follow-up. --- localization/assignCompartments.m | 383 +++++++++++++++++++ testing/function_tests/tAssignCompartments.m | 64 ++++ 2 files changed, 447 insertions(+) create mode 100644 localization/assignCompartments.m create mode 100644 testing/function_tests/tAssignCompartments.m diff --git a/localization/assignCompartments.m b/localization/assignCompartments.m new file mode 100644 index 00000000..37cb08de --- /dev/null +++ b/localization/assignCompartments.m @@ -0,0 +1,383 @@ +function [outModel, placement, addedTransports, exitFlag] = assignCompartments(model, GSS, reactionsToRelocate, varargin) +% assignCompartments Assign reactions to compartments by a functionality-constrained MILP. +% +% Deterministic alternative to predictLocalization: a single MILP places the requested reactions +% into compartments to agree with localization scores (GSS) WHILE keeping the model's objective +% (biomass) producible. Because functionality is a hard constraint, a reaction is placed against +% its own top score whenever the network needs it there (pathway coherence emerges from requiring +% flux, not from a heuristic). Mono-localization: each reaction is placed in exactly one +% compartment; a gene still spans compartments when it catalyses reactions placed in different ones. +% +% The model is merged to a single compartment first (like predictLocalization); the MILP then +% re-distributes the relocatable reactions across the compartments named in GSS, adding passive +% transports (defaultCompartment <-> c) where needed. +% +% Parameters +% ---------- +% model : struct +% a RAVEN model with an objective (model.c) and grRules. Multiple compartments are merged. +% GSS : struct +% gene scoring structure (genes, compartments, scores) as from parseScores. GSS.compartments +% are the target compartment labels and must include defaultCompartment. +% reactionsToRelocate : cell +% reaction ids to (re)place. Boundary reactions and the objective reaction are always pinned. +% +% Name-Value Arguments +% -------------------- +% defaultCompartment : char +% compartment that transports route through (usually cytosol); must be in GSS.compartments. +% transportCost : double (default 0.5) +% cost per added inter-compartment transport. +% multiCompartmentPenalty : double (default 0.5) +% cost per extra compartment a gene ends up in. +% minGrowth : double (default [] = 10%% of the unconstrained optimum) +% required objective flux. +% transportable : cell (default [] = all) +% metabolite ids (merged-model names) that may receive transports. Restricting it forces +% functionality-driven placement of reactions whose substrates are confined. +% bigM : double (default 1000) +% Big-M for flux gating. +% verbose : logical (default true) +% +% Returns +% ------- +% outModel : struct +% the compartmentalised model (objective still producible). +% placement : struct +% .rxns and .compartment (the assigned compartment id per relocated reaction). +% addedTransports : struct +% .mets and .compartment for the transports the MILP added. +% exitFlag : double +% 1 = optimal, -1 = infeasible/failed. +% +% See also +% -------- +% predictLocalization, parseScores, gapFillMILP + +p = parseRAVENargs(varargin, {'defaultCompartment',[]; 'transportCost',0.5; ... + 'multiCompartmentPenalty',0.5; 'minGrowth',[]; 'transportable',[]; 'bigM',1000; 'verbose',true}); +defaultCompartment = char(p.defaultCompartment); +transportCost = p.transportCost; +multiPen = p.multiCompartmentPenalty; +minGrowth = p.minGrowth; +M = p.bigM; +verbose = p.verbose; + +outModel = model; placement = struct('rxns',{{}},'compartment',{{}}); +addedTransports = struct('mets',{{}},'compartment',{{}}); exitFlag = -1; + +if all(model.c == 0) + error('RAVEN:badInput','model has no objective (set model.c).'); +end +comps = GSS.compartments(:); +if ~ismember(defaultCompartment, comps) + error('RAVEN:badInput','defaultCompartment ''%s'' not in GSS.compartments.', defaultCompartment); +end + +% Merge to a single compartment so every metabolite has one identity; relocations and transports +% are then expressed relative to defaultCompartment. +if numel(model.comps) > 1 + model = mergeCompartments(model, true, true); +end +nMet = numel(model.mets); +biomassIdx = find(model.c ~= 0, 1); + +% Determine min growth from the merged (single-compartment) model +if isempty(minGrowth) + sol = solveLP(model); + if isempty(sol.f) || sol.f <= 0 + error('RAVEN:badInput','merged model does not grow; pass minGrowth.'); + end + minGrowth = 0.1 * sol.f; +end + +% ---- scope: movable vs pinned ---- +isBoundary = (sum(model.S ~= 0, 1)' == 1); % one-metabolite reactions +relocate = ismember(model.rxns, reactionsToRelocate); +movableMask = relocate & ~isBoundary; +movableMask(biomassIdx) = false; +movIdx = find(movableMask); +pinIdx = find(~movableMask); +nMov = numel(movIdx); nPin = numel(pinIdx); nC = numel(comps); +defC = find(strcmp(comps, defaultCompartment)); + +% Genes in scope: those on movable reactions that have a row in GSS +[gInGSS, gssRow] = ismember(model.genes, GSS.genes); +geneOnMov = any(model.rxnGeneMat(movIdx, :) ~= 0, 1)'; +scopeMask = gInGSS & geneOnMov; +geneIdx = find(scopeMask); +nGene = numel(geneIdx); +score = zeros(nGene, nC); % gene x compartment score, columns aligned to comps +for gi = 1:nGene + score(gi, :) = GSS.scores(gssRow(geneIdx(gi)), :); +end + +% Transportable metabolites (touched by movable reactions) +movMetMask = any(model.S(:, movIdx) ~= 0, 2); +if isnumeric(p.transportable) && isempty(p.transportable) + transpMet = find(movMetMask); % default ([]): all movable metabolites transportable +else + transpMet = find(movMetMask & ismember(model.mets, p.transportable)); % a (possibly empty) set +end +% transport variables only for non-default compartments +trPairs = []; % [metRow, compCol] +for ci = 1:nC + if ci == defC; continue; end + trPairs = [trPairs; [transpMet, repmat(ci, numel(transpMet), 1)]]; %#ok +end +nTr = size(trPairs, 1); + +if verbose + fprintf('assignCompartments: %d movable, %d pinned reactions, %d genes, %d compartments, %d transports.\n', ... + nMov, nPin, nGene, nC, nTr); +end + +% ---- variable layout (columns) ---- +% fmove[mov,c] : nMov*nC continuous +% fpin[pin] : nPin continuous +% ftr[tr] : nTr continuous +% x[mov,c] : nMov*nC binary +% y[gene,c] : nGene*nC binary +% t[tr] : nTr binary +oFmove = 0; nFmove = nMov*nC; +oFpin = oFmove+nFmove; nFpin = nPin; +oFtr = oFpin+nFpin; % nTr +oX = oFtr+nTr; nX = nMov*nC; +oY = oX+nX; nY = nGene*nC; +oT = oY+nY; % nTr +nVar = oT+nTr; +fmoveCol = @(mi,ci) oFmove + (mi-1)*nC + ci; +xCol = @(mi,ci) oX + (mi-1)*nC + ci; +yCol = @(gi,ci) oY + (gi-1)*nC + ci; + +% ---- bounds ---- +lb = zeros(nVar,1); ub = zeros(nVar,1); +for k = 1:nMov + r = movIdx(k); + for ci = 1:nC + lb(fmoveCol(k,ci)) = min(model.lb(r),0); ub(fmoveCol(k,ci)) = max(model.ub(r),0); + end +end +lb(oFpin+(1:nPin)) = model.lb(pinIdx); ub(oFpin+(1:nPin)) = model.ub(pinIdx); +lb(oFtr+(1:nTr)) = -M; ub(oFtr+(1:nTr)) = M; +lb(oX+(1:nX)) = 0; ub(oX+(1:nX)) = 1; +lb(oY+(1:nY)) = 0; ub(oY+(1:nY)) = 1; +lb(oT+(1:nTr)) = 0; ub(oT+(1:nTr)) = 1; + +% ---- node balance: S over (met, compartment) ---- +% node row index = (metRow-1)*nC + ci +nNode = nMet*nC; +nodeRow = @(mr,ci) (mr-1)*nC + ci; +ri = []; ci_ = []; vv = []; +% movable: metabolites placed in c +for k = 1:nMov + r = movIdx(k); + mrows = find(model.S(:,r) ~= 0); + for ci = 1:nC + for mm = mrows' + ri(end+1)=nodeRow(mm,ci); ci_(end+1)=fmoveCol(k,ci); vv(end+1)=model.S(mm,r); %#ok + end + end +end +% pinned: metabolites in defaultCompartment +for k = 1:nPin + r = pinIdx(k); + mrows = find(model.S(:,r) ~= 0); + for mm = mrows' + ri(end+1)=nodeRow(mm,defC); ci_(end+1)=oFpin+k; vv(end+1)=model.S(mm,r); %#ok + end +end +% transports: -1 at (met,default), +1 at (met,c) +for k = 1:nTr + mm = trPairs(k,1); cc = trPairs(k,2); + ri(end+1)=nodeRow(mm,defC); ci_(end+1)=oFtr+k; vv(end+1)=-1; %#ok + ri(end+1)=nodeRow(mm,cc); ci_(end+1)=oFtr+k; vv(end+1)=+1; %#ok +end +A_node = sparse(ri, ci_, vv, nNode, nVar); +keepNode = any(A_node ~= 0, 2); % drop empty nodes +A_node = A_node(keepNode, :); +nNodeKept = size(A_node,1); + +% ---- flux gating: lb*x <= fmove <= ub*x ---- +giR=[]; giC=[]; giV=[]; row=0; gUb=[]; gLb=[]; +for k=1:nMov + r=movIdx(k); + for ci=1:nC + row=row+1; % ub: fmove - ub*x <= 0 + giR(end+1)=row; giC(end+1)=fmoveCol(k,ci); giV(end+1)=1; %#ok + giR(end+1)=row; giC(end+1)=xCol(k,ci); giV(end+1)=-max(model.ub(r),0); %#ok + gUb(end+1)=0; %#ok + end +end +A_gub = sparse(giR,giC,giV,row,nVar); b_gub=zeros(row,1); +giR=[]; giC=[]; giV=[]; row=0; +for k=1:nMov + r=movIdx(k); + for ci=1:nC + row=row+1; % lb: fmove - lb*x >= 0 + giR(end+1)=row; giC(end+1)=fmoveCol(k,ci); giV(end+1)=1; %#ok + giR(end+1)=row; giC(end+1)=xCol(k,ci); giV(end+1)=-min(model.lb(r),0); %#ok + end +end +A_glb = sparse(giR,giC,giV,row,nVar); b_glb=zeros(row,1); + +% ---- transport gating: -M*t <= ftr <= M*t ---- +tR=[];tC=[];tV=[]; +for k=1:nTr + tR(end+1)=k; tC(end+1)=oFtr+k; tV(end+1)=1; tR(end+1)=k; tC(end+1)=oT+k; tV(end+1)=-M; %#ok +end +A_tub = sparse(tR,tC,tV,nTr,nVar); % ftr - M*t <= 0 +tR=[];tC=[];tV=[]; +for k=1:nTr + tR(end+1)=k; tC(end+1)=oFtr+k; tV(end+1)=1; tR(end+1)=k; tC(end+1)=oT+k; tV(end+1)=M; %#ok +end +A_tlb = sparse(tR,tC,tV,nTr,nVar); % ftr + M*t >= 0 + +% ---- placement: sum_c x[mov,c] = 1 ---- +pR=[];pC=[];pV=[]; +for k=1:nMov + for ci=1:nC; pR(end+1)=k; pC(end+1)=xCol(k,ci); pV(end+1)=1; end %#ok +end +A_place=sparse(pR,pC,pV,nMov,nVar); b_place=ones(nMov,1); + +% ---- gene coupling x[r,c] <= y[g,c]; gene-has y<=sum x; gene1 sum_c y>=1 ---- +cR=[];cC=[];cV=[];row=0; +for k=1:nMov + gOn=find(model.rxnGeneMat(movIdx(k),:)~=0); + for g=gOn + gi=find(geneIdx==g,1); if isempty(gi); continue; end + for ci=1:nC + row=row+1; cR(end+1)=row; cC(end+1)=xCol(k,ci); cV(end+1)=1; + cR(end+1)=row; cC(end+1)=yCol(gi,ci); cV(end+1)=-1; %#ok + end + end +end +A_couple=sparse(cR,cC,cV,row,nVar); b_couple=zeros(row,1); +% gene-has: y[g,c] - sum_{r of g} x[r,c] <= 0 +hR=[];hC=[];hV=[];row=0; +for gi=1:nGene + g=geneIdx(gi); + rOfG=find(model.rxnGeneMat(movIdx,g)~=0)'; % indices into movIdx + for ci=1:nC + row=row+1; hR(end+1)=row; hC(end+1)=yCol(gi,ci); hV(end+1)=1; + for k=rOfG; hR(end+1)=row; hC(end+1)=xCol(k,ci); hV(end+1)=-1; end %#ok + end +end +A_has=sparse(hR,hC,hV,row,nVar); b_has=zeros(row,1); +% gene1: sum_c y[g,c] >= 1 +g1R=[];g1C=[];g1V=[]; +for gi=1:nGene + for ci=1:nC; g1R(end+1)=gi; g1C(end+1)=yCol(gi,ci); g1V(end+1)=1; end %#ok +end +A_gene1=sparse(g1R,g1C,g1V,nGene,nVar); b_gene1=ones(nGene,1); + +% ---- growth: fpin[biomass] >= minGrowth ---- +pinBio = find(pinIdx==biomassIdx,1); +A_grow = sparse(1, oFpin+pinBio, 1, 1, nVar); b_grow = minGrowth; + +% ---- objective: max sum score*y - multiPen*sum y - transportCost*sum t ---- +c = zeros(nVar,1); +for gi=1:nGene + for ci=1:nC; c(yCol(gi,ci)) = c(yCol(gi,ci)) + score(gi,ci) - multiPen; end +end +if isscalar(transportCost); tcost = repmat(transportCost,nTr,1); else; tcost = transportCost(trPairs(:,1)); end +c(oT+(1:nTr)) = c(oT+(1:nTr)) - tcost(:); + +% ---- assemble ---- +prob.A = [A_node; A_gub; A_glb; A_tub; A_tlb; A_place; A_couple; A_has; A_gene1; A_grow]; +prob.a = prob.A; +prob.b = [zeros(nNodeKept,1); b_gub; b_glb; zeros(nTr,1); zeros(nTr,1); b_place; b_couple; b_has; b_gene1; b_grow]; +prob.csense = [repmat('E',1,nNodeKept), repmat('L',1,numel(b_gub)), repmat('G',1,numel(b_glb)), ... + repmat('L',1,nTr), repmat('G',1,nTr), repmat('E',1,nMov), ... + repmat('L',1,numel(b_couple)), repmat('L',1,numel(b_has)), repmat('G',1,nGene), 'G']; +prob.c = -c; % optimizeProb minimises; we maximise c +prob.osense = 1; +prob.lb = lb; prob.ub = ub; +prob.vartype = repmat('C',1,nVar); +prob.vartype(oX+(1:nX)) = 'B'; +prob.vartype(oY+(1:nY)) = 'B'; +prob.vartype(oT+(1:nTr)) = 'B'; + +params.intTol = 1e-9; params.TimeLimit = 1000; +sol = optimizeProb(prob, params, verbose); +if ~checkSolution(sol) + if verbose; fprintf('assignCompartments: MILP infeasible or failed.\n'); end + return; +end + +% ---- extract ---- +xval = sol.full(oX+(1:nX)); +placeRxns = {}; placeComp = {}; +for k=1:nMov + [~,ci] = max(xval((k-1)*nC + (1:nC))); + placeRxns{end+1} = model.rxns{movIdx(k)}; placeComp{end+1} = comps{ci}; %#ok +end +placement.rxns = placeRxns(:); placement.compartment = placeComp(:); +tval = sol.full(oT+(1:nTr)); +sel = find(tval > 0.5); +addedTransports.mets = model.mets(trPairs(sel,1)); +addedTransports.compartment = comps(trPairs(sel,2)); +exitFlag = 1; + +outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment); +end + +% ------------------------------------------------------------------- apply +function outModel = i_applyAssignment(model, placement, addedTransports, comps, defaultCompartment) +% Build the compartmentalised model: relabel each placed reaction's metabolites into its +% compartment (creating per-compartment metabolite copies) and add the requested transports. +outModel = model; +% ensure all target compartments exist +for ci = 1:numel(comps) + if ~ismember(comps{ci}, outModel.comps) + outModel.comps{end+1,1} = comps{ci}; + if isfield(outModel,'compNames'); outModel.compNames{end+1,1} = comps{ci}; end + end +end +% metabolite (baseName, compartment) -> index, seeded with existing +metKey = strcat(model.metNames, '###', model.comps(model.metComps)); +for k = 1:numel(placement.rxns) + r = find(strcmp(outModel.rxns, placement.rxns{k}), 1); + comp = placement.compartment{k}; + mrows = find(model.S(:, find(strcmp(model.rxns,placement.rxns{k}),1)) ~= 0); %#ok<*FNDSB> + for mm = mrows' + [outModel, newMet] = i_metInComp(outModel, model, mm, comp); + outModel.S(newMet, r) = model.S(mm, find(strcmp(model.rxns,placement.rxns{k}),1)); + if ~isequal(newMet, mm); outModel.S(mm, r) = 0; end + end +end +for k = 1:numel(addedTransports.mets) + outModel = i_addTransport(outModel, model, addedTransports.mets{k}, addedTransports.compartment{k}, defaultCompartment); +end +end + +function [outModel, idx] = i_metInComp(outModel, model, srcMetRow, comp) +ci = find(strcmp(outModel.comps, comp), 1); +name = model.metNames{srcMetRow}; +cand = find(strcmp(outModel.metNames, name) & outModel.metComps == ci, 1); +if ~isempty(cand); idx = cand; return; end % reuse existing per-compartment metabolite +newId = [model.mets{srcMetRow} '_' comp]; +outModel.mets{end+1,1} = newId; +outModel.metNames{end+1,1} = name; +outModel.metComps(end+1,1) = ci; +outModel.S(end+1, :) = 0; +if isfield(outModel,'b'); outModel.b(end+1,1) = 0; end +idx = numel(outModel.mets); +end + +function outModel = i_addTransport(outModel, model, metId, comp, defaultCompartment) +srcRow = find(strcmp(model.mets, metId), 1); +if isempty(srcRow); return; end +[outModel, dRow] = i_metInComp(outModel, model, srcRow, defaultCompartment); +[outModel, cRow] = i_metInComp(outModel, model, srcRow, comp); +trId = ['tr_' model.metNames{srcRow} '_' comp]; +if any(strcmp(outModel.rxns, trId)); return; end +outModel.rxns{end+1,1} = trId; +outModel.S(:, end+1) = 0; outModel.S(dRow, end) = -1; outModel.S(cRow, end) = 1; +outModel.lb(end+1,1) = -1000; outModel.ub(end+1,1) = 1000; +if isfield(outModel,'rev'); outModel.rev(end+1,1) = 1; end +if isfield(outModel,'c'); outModel.c(end+1,1) = 0; end +if isfield(outModel,'rxnNames'); outModel.rxnNames{end+1,1} = trId; end +if isfield(outModel,'grRules'); outModel.grRules{end+1,1} = ''; end +if isfield(outModel,'rxnGeneMat'); outModel.rxnGeneMat(end+1,:) = 0; end +end diff --git a/testing/function_tests/tAssignCompartments.m b/testing/function_tests/tAssignCompartments.m new file mode 100644 index 00000000..756352c2 --- /dev/null +++ b/testing/function_tests/tAssignCompartments.m @@ -0,0 +1,64 @@ +classdef tAssignCompartments < RavenTestCase +% tAssignCompartments Tests for the functionality-constrained compartment-assignment MILP. +% +% Uses a tiny single-compartment toy: EX_A (uptake) -> r1: A->B -> bio: B-> (objective), with +% gene g1 (on r1) scored higher for a hypothetical compartment 'm' than for the cytosol 'c'. +% The decisive checks are that the functionality constraint can override the localization score +% and that the returned model stays able to grow. Requires a MILP solver. + + methods (Test) + + function functionalityOverridesScore(testCase) + testCase.assumeMILPSolver(); + model = tAssignCompartments.toy(); + GSS = tAssignCompartments.gss(); + % B is non-transportable, so placing r1 in 'm' would strand B where biomass (in 'c') + % cannot reach it -> r1 must stay in 'c', against g1's higher 'm' score. + evalc(['[oModel, place, ~, eFlag] = assignCompartments(model, GSS, {''r1''}, ' ... + '''defaultCompartment'', ''c'', ''transportable'', {}, ''verbose'', false);']); + testCase.verifyEqual(eFlag, 1); + testCase.verifyEqual(place.compartment{1}, 'c'); + sol = solveLP(oModel); + testCase.verifyGreaterThan(sol.f, 1e-6); % applied model still grows + end + + function scoreRespectedWhenFunctional(testCase) + testCase.assumeMILPSolver(); + model = tAssignCompartments.toy(); + GSS = tAssignCompartments.gss(); + % Transports allowed and cheap -> r1 can follow its score to 'm' and stay functional. + evalc(['[oModel, place, trans, eFlag] = assignCompartments(model, GSS, {''r1''}, ' ... + '''defaultCompartment'', ''c'', ''transportCost'', 0.1, ''verbose'', false);']); + testCase.verifyEqual(eFlag, 1); + testCase.verifyEqual(place.compartment{1}, 'm'); + testCase.verifyGreaterThan(numel(trans.mets), 0); % transports were added + sol = solveLP(oModel); + testCase.verifyGreaterThan(sol.f, 1e-6); + end + + function infeasibleGrowthFloorReported(testCase) + testCase.assumeMILPSolver(); + model = tAssignCompartments.toy(); + GSS = tAssignCompartments.gss(); + evalc(['[~, ~, ~, eFlag] = assignCompartments(model, GSS, {''r1''}, ' ... + '''defaultCompartment'', ''c'', ''minGrowth'', 1e6, ''verbose'', false);']); + testCase.verifyEqual(eFlag, -1); % growth floor unreachable + end + + end + + methods (Static) + function model = toy() + model.id='toy'; model.comps={'c'}; model.compNames={'cytoplasm'}; + model.mets={'A_c';'B_c'}; model.metNames={'A';'B'}; model.metComps=[1;1]; + model.S=sparse([-1 -1 0; 0 1 -1]); % [A;B] x [EX_A r1 bio] + model.rxns={'EX_A';'r1';'bio'}; model.rxnNames=model.rxns; + model.lb=[-10;0;0]; model.ub=[1000;1000;1000]; model.rev=[1;0;0]; + model.c=[0;0;1]; model.b=zeros(2,1); + model.genes={'g1'}; model.grRules={'';'g1';''}; model.rxnGeneMat=sparse([0;1;0]); + end + function GSS = gss() + GSS.genes={'g1'}; GSS.compartments={'c';'m'}; GSS.scores=[0.4 0.9]; + end + end +end