demo_yeast_ccm_reconstruct_model

 DEMO_YEAST_CCM_RECONSTUCT_MODEL - Demo script for reconstruction of enzyme-balanced models

0001 % DEMO_YEAST_CCM_RECONSTUCT_MODEL - Demo script for reconstruction of enzyme-balanced models
0002 
0003 demo_dir = [fileparts(which(mfilename))];
0004 
0005 cd(demo_dir)
0006 
0007 echo on;
0008 clc
0009 %---------------------------------------------------------------------------------
0010 % DEMO: Economic flux analysis for yeast central metabolism
0011 %
0012 % In this script, we build an enzyme-balanced kinetic model in the following way:
0013 %
0014 % 1.  Determine flux distribution using principle of minimal fluxes
0015 % 2a. Choose thermodynamic forces
0016 % 2b. Choose economic potentials (with the principle of homogeneous costs)
0017 % 3.  Reconstruct an economic-kinetic model realising the flux distribution
0018 %---------------------------------------------------------------------------------
0019  
0020 % Press key to continue
0021  
0022 pause
0023 clc
0024 % --------------------------------------------------------------------
0025 % We load a network model of yeast central carbon metabolism.
0026  
0027 % The file contains the variables: network, network_CoHid, network_CoSplit, v_sign
0028  
0029 load('/data/yeast_ccm_network.mat');
0030  
0031 % Press key to continue
0032  
0033 pause
0034 clc
0035 % --------------------------------------------------------------------
0036 % Phase 1: Determine flux distribution (principle of minimal fluxes)
0037 %          for ATP production, with some predefined flux directions
0038 % ------------------------------------------------------------------
0039  
0040 % We create structs 'cba_options' and 'cba_constraints' with some default settings:
0041  
0042 [cba_options, cba_constraints] = cba_default_options(network);
0043  
0044 % Press key to continue
0045  
0046 pause
0047 clc
0048 % ------------------------------------------------------------------
0049 % We adjust the settings for flux constraints
0050  
0051 cba_constraints.v_sign = v_sign;
0052  
0053 % We define some reactions to be inactive
0054 % (strings like 'R00253' are reaction names in our network)
0055  
0056 cba_constraints.v_fix(label_names({'R00253'},network.actions)) = 0;
0057 cba_constraints.v_fix(label_names({'R00258'},network.actions)) = 0;
0058 cba_constraints.v_fix(label_names({'R00114'},network.actions)) = 0;
0059 cba_constraints.v_fix(label_names({'R00243'},network.actions)) = 0;
0060 cba_constraints.v_fix(label_names({'R00341'},network.actions)) = 0;
0061 cba_constraints.v_fix(label_names({'R00344'},network.actions)) = 0;
0062 cba_constraints.v_fix(label_names({'R00711'},network.actions)) = 0;
0063  
0064 % We fix one flux direction
0065  
0066 cba_constraints.v_sign(label_names({'R00342'},network.actions)) = 1;
0067  
0068 % We set an upper bound for one reaction
0069  
0070 cba_constraints.v_max(label_names({'Oxphos__NADH__irrev__ATP'},network.actions)) = 2;
0071  
0072 % Press key to continue
0073  
0074 pause
0075 clc
0076 % ------------------------------------------------------------------
0077 % Next, we adjust the settings for the metabolic objective
0078  
0079 % Our metabolic objective is the net production of ATP.
0080 % We declare this by a setting a vector zx with weights
0081 % for all external production rates
0082  
0083 zx = zeros(nm,1); 
0084 zx(label_names({'ATP'},network.metabolites)) = 1;
0085  
0086 % Using this vector, we define the marginal gain vectors (z_ext and zv)
0087  
0088 cba_constraints.z_int = 0 * cba_constraints.z_int;
0089 cba_constraints.z_ext = zx(find(network.external));
0090 cba_constraints.zv    = network.N' * zx;
0091  
0092 % Press key to continue
0093  
0094 pause
0095 clc
0096 % ------------------------------------------------------------------
0097 % Now we determine the stationary fluxes
0098  
0099 
0100 % 1. Flux Balance Analysis
0101  
0102 [v_fba,f_benefit] = fba(network,cba_constraints);
0103  
0104  
0105 % 2. Fix objective achieved in FBA and run Flux Minimisation (PMF)
0106  
0107 f_benefit = cba_constraints.zv'*v_fba; 
0108  
0109 v = pmf(network,cba_constraints,f_benefit,v_fba);
0110  
0111  
0112 % 3. Set very small fluxes to 0 and make the flux mode stationary
0113  
0114 v(abs(v) < 10^-5 *max(abs(v))) = 0;
0115  
0116 v = project_fluxes(network.N,find(network.external), v,[],sign(v),struct('method','euclidean'));
0117  
0118  
0119 % Press key to continue
0120  
0121 pause
0122 clc
0123 % ---------------------------------------------------------------
0124 % Just to be sure, we apply a function that would remove
0125 % thermodynamically infeasible cycles from our flux mode
0126  
0127 cba_constraints.ind_ignore = label_names({'Biomass_production'},network.actions);
0128  
0129 [v,C] = eba_make_feasible(v, network. N, 'loose', nan, cba_constraints.ind_ignore);
0130  
0131 % Press key to continue
0132  
0133 pause
0134 clc
0135 % --------------------------------------------------------------------
0136 % Phase 2: Choose chemical and economic potentials
0137  
0138 % --------------------------------------------------------------------
0139 % We adjust the settings for choosing the chemical potentials
0140  
0141 cba_constraints.dmu_limit     = 10;
0142 cba_constraints.dmu_limit_min = 2;
0143 cba_constraints.mu_min        = -20 * ones(size(cba_constraints.mu_min));
0144 cba_constraints.mu_max        =  20 * ones(size(cba_constraints.mu_min));
0145 cba_constraints.dmu_min       = -20 * ones(size(cba_constraints.dmu_min));
0146 cba_constraints.dmu_max       =  20 * ones(size(cba_constraints.dmu_min));
0147 cba_constraints.rho           = 100;
0148   
0149 % Press key to continue
0150  
0151 pause
0152 clc
0153 % --------------------------------------------------------------------
0154 % We choose the chemical potentials
0155  
0156 [mu, success_flag] = sample_feasible_mu(network.N,find(network.external),v,cba_constraints,cba_options,'sample',1);
0157   
0158 % Press key to continue
0159  
0160 pause
0161 clc
0162 % --------------------------------------------------------------------
0163 % We choose the economic potentials
0164  
0165 [w, delta_w, y, zx] = cba_homogeneous_cost(network, v, cba_constraints);
0166  
0167 % Press key to continue
0168  
0169 pause
0170 clc
0171 % --------------------------------------------------------------------
0172 % Phase 3: Reconstruct an enzyme-balanced model
0173  
0174 % We set metabolite levels. For the sake of this example, we simply use random values
0175   
0176 c = 1+5*rand(size(mu));
0177  
0178 % Now we set some options for the model reconstruction ..
0179  
0180 cba_options.check_curvatures = 0; 
0181 cba_constraints.Q_ext        = [];
0182 cba_constraints.hu           = ones(size(network.actions)); 
0183  
0184 % .. and reconstruct the model
0185  
0186 [network_new, res, cba_constraints_new] = cba_reconstruct_model(network, v, mu, cba_constraints, cba_options, y, w, c);
0187   
0188 % Press key to continue
0189  
0190 pause
0191 clc
0192 % ------------------------------------------------------------------
0193 % Finally, we plot some of the reconstructed quantities
0194  
0195 % 1. Network and external metabolites
0196  
0197 figure(1); clf; 
0198 netgraph_concentrations(network_CoSplit,network.external,[],1,struct('actprintnames',1));
0199   
0200 % Press key to continue
0201 pause
0202  
0203 % 2. Metabolic fluxes
0204  
0205 figure(2); clf; 
0206 netgraph_concentrations(network_CoSplit,-network.external,v,1,struct('actstyle','none','arrowsize',0.03));
0207   
0208 % Press key to continue
0209 pause
0210  
0211 % 3. Predefined flux signs
0212  
0213 figure(3); clf; 
0214 netgraph_concentrations(network_CoHid,[],v_sign,1,struct('actstyle','fixed','arrowsize',0.03,'actprintnames',1));
0215   
0216 % Press key to continue
0217 pause
0218  
0219 % 4. The local ATP production (contribution of reactions to the metabolic gain)
0220  
0221 figure(4); clf; 
0222 netgraph_concentrations(network_CoSplit,[],cba_constraints.zv .* v,1,struct('arrowstyle','none'));
0223   
0224 % Press key to continue
0225 pause
0226  
0227 % 5. The chemical potentials
0228  
0229 figure(5); clf;
0230 netgraph_concentrations(network_CoHid,mu,[-network.N'*mu].*[v~=0],1,struct('actstyle','none','arrowsize',0.01));
0231   
0232 % Press key to continue
0233 pause
0234  
0235 % 6. The economic potentials
0236  
0237 figure(5); clf;
0238 netgraph_concentrations(network_CoHid,w,[-network.N'*w].*[v~=0],1,struct('actstyle','none','arrowsize',0.01));
0239  
0240 % Press key to continue
0241 pause
0242 clc
0243 % That was it - we built a nice kinetic model!
0244  
0245 % All variables are still in your workspace.
0246  
0247 % Enjoy working with elasticity sampling!
0248  
0249 % Press key to finish
0250 pause
0251 return