% Add your path to cplex to the matlab path  
% Windows e.g.: addpath(genpath('C:\Program Files\IBM\ILOG\CPLEX_Studio1271'))
% MAC e.g: addpath(genpath('/Users/USERNAME/Applications/IBM/ILOG/CPLEX_Studio127')) 
% Add your path to mattfa toolbox in the matlab path  
% addpath('C:\path\to\mattfa')

% Change the solver 
changeCobraSolver('cplex_direct')

%% Load the model
tmp = load('smallEcoli.mat');
model = tmp.smallEcoli;
clear tmp

%% Load the thermodynamics database
tmp = load('thermo_data.mat');
ReactionDB = tmp.DB_AlbertyUpdate;
clear tmp


%% Perfom an FBA on the model (optimize for growth)
solFBA = optimizeCbModel(model);

%% Prepare the model for TFA
prepped_model = prepModelforTFA(model, ReactionDB, model.CompartmentData);

%% Convert to TFA
min_obj = 0.0
tmp = convToTFA(prepped_model, ReactionDB, [], 'DGo', [], min_obj);

% Add net flux variables, which are equal to forwards flux - backwards flux
% The netflux variables work similar to the FBA fluxes
% NF_rxn = F_rxn - B_rxn
this_tmodel = addNetFluxVariables(tmp);

%% Perfom TFA on the model (optimise for growth)
soltFA = solveTFAmodelCplex(this_tmodel);


%% We add some generic data for concentrations
metNames = {'adp_c', 'atp_c'};
C_lb = [1e-06, 1e-03]';
C_ub = [7e-04, 5e-02]';
LC_varNames = {'LC_adp_c', 'LC_atp_c'};
% find the indices of these variables in the variable names of the tfa

id_LC_varNames = find_cell(LC_varNames, this_tmodel.varNames);
% Set to the model these log-concentration values
this_tmodel.var_lb(id_LC_varNames) = log(C_lb);
this_tmodel.var_ub(id_LC_varNames) = log(C_ub)

%% Optimise for growth with concentration data 
soltFA_w_concentrations = solveTFAmodelCplex(this_tmodel)