Model analysis¶

Once a model ensemble has been generated by GRASP, we can easily:

simulate the models;

do metabolic control analysis;

do stability analysis.

Simulating the model ensemble¶

To simulate the model you can use the script simulate_model.m in the examples folder.

It is worth keeping in mind that GRASP uses scaled enzyme and metabolite concentrations. This is important when simulating perturbations.

While for enzymes, the concentrations are always scaled concentrations and any perturbations applied will be relative to the reference state, for metabolites perturbations can be defined either relative to the reference state or in terms of absolute concentrations, by setting the variable metsAbsOrRel to 'rel' or 'abs', respectively.

As an example, if the initial concentration for a given enzyme is set to 2, this means its concentration will be twice its concentration in the reference state. The same is true for metabolites if metsAbsOrRel is set to 'rel'.

If you want to define the initial concentration for a given metabolite in terms of absolute concentrations, metsAbsOrRel must be set to 'abs' and the initial concentration must be specified in mol/L.

Other important variables to define are:

finalTime: the time of the simulation given in hours;

interruptTime: when the simulation will be aborted if it is still running, this time is specified in seconds. The idea is to skip simulations that take a long time to run (probably something is wrong with them anyways);

Example:

% Add all GRASP functions to the Matlab path
addpath(fullfile('..', 'matlab_code', 'analysisFxns'), ...
        fullfile('..', 'matlab_code', 'ensembleFxns'), ...
        fullfile('..', 'matlab_code', 'patternFxns'));


% Define the model ID (file name of the model ensemble)
modelID = 'toy_model';

% Define the path to the input/output files
outputFolder = fullfile(, '..', 'io', 'output');

% How many models in the ensemble you want to simulate
numModels = 5;

% How many cores you want to use to run the simulation
numCores = 2;

% Define how many seconds until ODE solver is interrupted. The idea is to
% skip models that take ages to simulate.
interruptTime = 40;

% Load model ensemble previously generated
load(fullfile(outputFolder, [modelID, '.mat']))

% Define whether the initial condition for metabolites is a relative or an
% absolute concentration by setting metsAbsOrRel to either 'rel' or 'abs',
% respectively
metsAbsOrRel = 'rel';

% Change initial conditions here if you want, format: {rxn/met ID, initial value}
enzymesIC = {{'r1', 1}, {'r10', 1}};       % Always relative concentrations for enzymes
metsIC = {{'m5', 0.6}, {'m11', 1.2}};      % When setting metsAbsOrRel = 'abs', absolute concentrations must be given in mol/L

% Specifiy the time of simulation (in hours)
finalTime = 1;

% Run simulation
simulationRes = simulateEnsemble(ensemble, finalTime, enzymesIC, metsIC, metsAbsOrRel, interruptTime, numModels, numCores);


% Save the results
save(fullfile(outputFolder, ['simulation_', modelID, '.mat']), 'simulationRes')
write(cell2table(ensemble.mets(ensemble.metsActive)), fullfile(outputFolder, [modelID, '_metsActive.dat']));
write(cell2table(ensemble.rxns(ensemble.activeRxns)), fullfile(outputFolder, [modelID, '_rxnsActive.dat']));

The resulting Matlab structure contains the following fields:

t: time points in each model simulation;

conc: metabolite concentrations (relative) for each time point and model simulation;

flux: reaction fluxes (absolute) for each time point and model simulation.

Metabolite concentrations are always returned as relative values while reaction fluxes are always returned as absolute values.

To visualize the simulation results you can use the jupyter notebook visualize_simulations.ipynb in the visualization folder. More details about this notebook can be found in the Visualizing results section.

Metabolic Control Analysis (MCA)¶

To do Metabolic Control Analysis (MCA) on the model ensemble you can use the MCA_analysis.m script in the examples folder.

This will perform MCA on the model ensemble and give the average flux and concentration control coefficients over all models in the ensemble.

If the variable saveMCAMatrices is set to 1 it wil also return the control coefficients for each individual model. This can be useful to do further analysis. However, due to memory requirements it can make the calculation very slow for larger models (> 40 reactions).

The resulting Matlab structure contains the following fields:

xControlAvg: average concentration control coefficient for each metabolite over the whole model ensemble;

vControlAvg: average flux control coefficient for each reaction over the whole model ensemble;

xcounter: number of models in the average concentration control coefficient calculation;

vcounter: number of models in the average flux control coefficient calculation;

xControl: concentration control coefficient matrix for each model;

vControl: flux control coefficient matrix for each model;

E_x_nor: normalized elasticity matrix for each model.

Example:

% Add all GRASP functions to Matlab's path
addpath(fullfile('..', 'matlab_code', 'patternFxns'), ...
        fullfile('..', 'matlab_code', 'ensembleFxns'));

% Whether or not to save the MCA results for all models and not just mean values
saveMCAMatrices = 1;

% Define the model ID (file name of the model ensemble)
modelID = 'toy_model';

% Define the path to the input/output files
outputFolder = fullfile(, '..', 'io', 'output');

% Load previously generated model ensemble
load(fullfile(outputFolder, [modelID, '.mat']))

% Run MCA analysis
mcaResults = controlAnalysis(ensemble, saveMCAMatrices);

% Save MCA results
save(fullfile(outputFolder, ['MCA_', modelID, '.mat']), 'mcaResults');
write(cell2table(ensemble.rxns(ensemble.activeRxns)), fullfile(outputFolder, [modelID, '_rxnsActive.dat']));
write(cell2table(ensemble.mets(ensemble.metsActive)), fullfile(outputFolder, [modelID, '_metsActive.dat']));
write(cell2table(mcaResults.enzNames), fullfile(outputFolder, [modelID, '_enzNames.dat']));

% Plot MCA results - optional

% Optional, Define ranges for displaying the MCA results``:  {1st category, range; 2nd category, range}
% For example, categories = {'Glycolysis',[1,20]; 'Pentose Phosphate Pathway',[21,30];'Others', [31,37]};
categories = {'all', [1, 5]};

plotControlAnalysis(mcaResults, ensemble, categories);

If you have promiscuous enzymes in your model, you should do response analysis instead of simple control analysis.

This is because an increase of the promiscuous enzyme concentration doesn’t necessarily lead to an equal increase in the flux of the reactions it catalyzes.

Response and control coefficients are the same when enzymes are independent and an increase in enzyme concentration leads to a proportional increase in the reaction flux.

Response coefficients are calculated as:

\[C_E^J = C_v^J\Pi\]

where \(C_v^J\) is the flux control coefficient matrix and \(\Pi\) is the parameter elasticity matrix.

The implementation is based on

Computational Models of Metabolism: Stability and Regulation in Metabolic Networks, section on MCA;

Enzyme-enzyme interactions and control analysis: 1. The case of non-additivity: monomer-oligomer associations;

Enzyme-enzyme interactions and control analysis: 2. The case of non-independence: heterologous associations.

You can do response analysis by using the functions controlAndResponseAnalysis instead of controlAnalysis and plotControlAndResponseAnalysis instead of plotControlAnalysis.

To visualize the results you can use the jupyter notebook visualize_mca.ipynb in the visualization folder. For more details see the Visualizing results section.

Stability analysis¶

To do stability analysis you can use the stability_analysis.m script in the examples folder.

This basically calculates the jacobian of every model in the ensemble and checks if the real part of its eigenvalues is higher than the given threshold, eigThreshold.

If the specified threshold is the same as the one specified when building the model ensemble, all models will be considered stable, since only stable models are returned by GRASP when building the model ensemble.

The returned Matlab structure has the following fields:

posEig: positive eigenvalues for unstable models;

unstableModels: list of unstable models.

Example:

% Add all GRASP functions to Matlab's path
addpath(fullfile('..', 'matlab_code', 'patternFxns'), ...
        fullfile('..', 'matlab_code', 'ensembleFxns'));

% Whether or not to save the MCA results for all models and not just mean values
saveMCAMatrices = 1;

% Define the model ID (file name of the model ensemble)
modelID = 'toy_model';

% Define the path to the input/output files
outputFolder = fullfile(, '..', 'io', 'output');

% threshold of the jacobian's eigenvalues
eigThreshold = 10^-5;

% Load the model ensemble generated by GRASP
load(fullfile(outputFolder, [modelID, '.mat']));

% Run stability analysis
stabilityRes = ensembleStabilityTest(ensemble, eigThreshold);

% Save the results
save(fullfile(outputFolder, ['stability_', modelID, '.mat']), 'stabilityRes');