Pscmtoolbox

• Bmrcalculations
• Hostmicrobeinteraction
• Setconstraints
  • Inputdata
  • Organweight
  • Diets
• Io

LeakTestRecon

changeCobraSolver(‘tomlab_cplex’,’lp’); modelClosed = modelConsistent;

Test4HumanFctExtv5(model, test, optionSinks)

function [TestSolution,TestSolutionName] = Test4HumanFct(model,test) This functions test for the ~460 human functions - I removed duplicates

INPUT model model structure (Recon1, with desired in silico

condition)

test possible statements: Recon1, IECori, IEC, all (default)

(choose IECori if you intend to test the IEC model OR a model that contains lumen (‘u’) as compartment otw choose IEC); all check for Recon1 and IEC

option if true = set sink reactions to 0 (default, leave unchanged).

Note that all lb’s of exchanges and demands will be set to 0

OUTPUT TestSolution array containing the optimal value for the different

tests

TestSolutionName array containing the names for the different tests

Ines Thiele, 09/05/09 MKA, 03/04/12 some of the reaction names have changed in newer versions of Recon1 and Recon2. Comment setup til line 146 if using an old version. MKA, 24/05/12 finds correct EX_reactions and changes these to zero IT, 07/20/12 added tests for sIEC model IT 2017 more tests added

adjustWholeBodyRxnCoeff(model, listOrgan, listCoeff)

[model] = adjustWholeBodyRxnCoeff(model, listOrgan, listCoeff)

This function adjusts the coefficients of the whole-body biomass maintenance (WBM) reaction. The WBM reaction contains each organ present in the whole-body metabolic reconstructions. For each organ, the stoichiometric coefficients represent the fractional weight contribution of the respective organ to the whole body weight. These coefficients can be updated to reflect individual specific body contributions. E.g., in obese individuals the ratio of muscle and adipose tissue is different than in a normal BMI individual. Hence, they can be updated with this function.

INPUT model whole-body metabolic model listOrgan List of organs, whose stoichiometric coefficient should be

updated

listCoeff List of coefficients that replace current ones in the WBM

reaction (order must match the order of organs in ListOrgan)

OUTPUT model whole-body metabolic model with adjusted stoichiometric

coefficients in the whole-body metabolic model

Ines Thiele, 2012 - 2020

annotateHH

annotate Harvey and Harvetta

annotateModel(model, annotateRxns, annotateMets, modelID, modelName, modelAnnotation)

This function annotates a model with VMH reaction and metabolite identifiers.

function model = annotateModel(model, annotateRxns,annotateMets)

INPUT model Model structure annotateRxns default: 1 annotateMets default: 1

Optional: modelID ID of model modelName model Name

OUTPUT model Updated model structure

Ines Thiele October 2019

annotateOrganAtlas

annotate OrganAtlas

calcOrganFract(model, IndividualParameters)

This function extrapolates the organ weight fractions based on polynomials given in http://www.ams.sunysb.edu/~hahn/psfile/pap_obesity.pdf (PMID 19267313), Table 3. Those ones that are not given are assumed to remain constant with weight using the fractions from the reference man and reference woman.

[organs,OrganWeight,OrganWeightFract,IndividualParameters] = calcOrganFract(model, IndividualParameters)

INPUT model Model structure IndividualParameters Structure of individual parameters (sex, weight,

blood volume)

OUTPUT organs List of organs OrganWeight List of organ weights (same order as organs) OrganWeightFract List of organ weight fractions (same order as organs) IndividualParameters Updated structure of individual parameters

Ines Thiele Nov 2017

calculateBMR(sex, weight, height, age)

This function calculates the basal metabolic rate using the phenomenological model proposed by Harris-Benedict and derivations thereof (see below. Please also refer to the corresponding wikipedia entry: https://en.wikipedia.org/wiki/Harris%E2%80%93Benedict_equation for more details

function BMR = calculateBMR(sex, weight, height, age)

INPUT sex ‘male’ or ‘female’ weight in kg height in cm age in years

OUTPUT BMR array with 3 numbers calculated based on

1. original Harris Benedict equations [1],[2]

2. Harris Benedict equations revised by Roza and Shizgal in 1984.[3]

3. The Harris Benedict equations revised by Mifflin and St Jeor in 1990:[4]

References: [1] Harris JA, Benedict FG (1918). “A Biometric Study of Human Basal Metabolism”. Proceedings of the National Academy of Sciences of the United States of America. 4 (12): 370?3. doi:10.1073/pnas.4.12.370. PMC 1091498?Freely accessible. PMID 16576330. [2] A Biometric Study of Basal Metabolism in Man. J. Arthur Harris and Francis G. Benedict. Washington, DC: Carnegie Institution, 1919. [3] Roza AM, Shizgal HM (1984). “The Harris Benedict equation reevaluated: resting energy requirements and the body cell mass”. The American Journal of Clinical Nutrition. 40 (1): 168?82. PMID 6741850. [4] Mifflin MD, St Jeor ST, Hill LA, Scott BJ, Daugherty SA, Koh YO (1990). “A new predictive equation for resting energy expenditure in healthy individuals”. The American Journal of Clinical Nutrition. 51 (2): 241?7. PMID 2305711.

Ines Thiele 01/2018

checkIEM_WBM

This function performs the inborn error of metabolism simulations by deleting (or reducing) the flux through reaction(s) affected by a gene defect and optimized the flux through a defined set of biomarker reactions.

function [IEMSol] = checkIEM_WBM(model,IEMRxns, BiomarkerRxns,minRxnsFluxHealthy, reverseDirObj, fractionKO,minBiomarker,fixIEMlb)

INPUT model whole-body metabolic reconstruction or Recon3Dmodel IEMRxns Reaction(s) affected by the inborn error of

metabolism

minRxnsFluxHealthy min flux value(s) through the IEMRxns reverseDirObj the function maximizes the objective flux by

default. If set to 1, the function also checks the minimization problem.

fractionKO By default, a complete knowckout of BiomarkerRxnsthe IEM

reactions is computed but it is possible to set a fraction (default = 1 for 100% knockout)

minBiomarker Minimization through biomarker reaction (default = 0) fixIEMlb fix IEM to lb = ub

=(1-fractionKO)*solution.v(find(model.c)) (default = 0, i.e., lb =0, while ub = (1-fractionKO)*solution.v(find(model.c))

LPSolver Define LPSolver (‘ILOGcomplex’ - default;

‘tomlab_cplex’)

OUTPUT IEMSol Predicted biomarker fluxes and comparison with the

reported biomarkers

USAGE: Exampe of preparation of a set of inputs to checkIEM_WBM

R = {‘_2OXOADPTm’;’_2AMADPTm’;’_r0879’}; RxnsAll2 = ‘’; for i = 1: length(R)

RxnsAll = model.rxns(find(~cellfun(@isempty,strfind(model.rxns,R{i})))); RxnsAll2 =[RxnsAll2;RxnsAll];

end IEMRxns = unique(RxnsAll2); RxnMic = model.rxns(find(~cellfun(@isempty,strfind(model.rxns,’Micro_’)))) ; IEMRxns = setdiff(IEMRxns,RxnMic);

if ~strcmp(modelName,’Recon3D’)

% add demand reactions to blood compartment for those biomarkers reported for blood % biomarker based on https://www.omim.org/entry/204750 model = addDemandReaction(model, ‘L2aadp[bc]’);

BiomarkerRxns = {

‘DM_L2aadp[bc]’ ‘Increased (blood)’ ‘EX_2oxoadp[u]’ ‘Increased (urine)’ ‘EX_adpoh[u]’ ‘Increased (urine)’ };

else

BiomarkerRxns = {

‘EX_2oxoadp[u]’ ‘Increased (urine)’ ‘EX_adpoh[u]’ ‘Increased (urine)’ };

end

Then call the checkIEM_WBM function

[IEMSol_2OAA] = checkIEM_WBM(model,IEMRxns, BiomarkerRxns,minRxnsFluxHealthy);

Example IEMSol returned from the above

{‘IEM Rxns All obj - Healthy’} {‘65403.5393’ } {0×0 double } {‘IEM Rxns All obj - Disease’} {‘65403.5393’ } {0×0 double } {‘WB obj - Healthy’ } {‘NA’ } {0×0 double } {‘WB obj - Disease’ } {‘1’ } {‘1’ } {‘Healthy:DM_L2aadp[bc]’ } {‘-5.5324e-08’} {‘Disease - Reported:Increas…’} {‘Disease:DM_L2aadp[bc]’ } {‘40.2418’ } {‘Disease - Reported:Increas…’} {‘Healthy:EX_2oxoadp[u]’ } {‘-2.3603e-10’} {‘Disease - Reported:Increas…’} {‘Disease:EX_2oxoadp[u]’ } {‘3.7368’ } {‘Disease - Reported:Increas…’} {‘Healthy:EX_adpoh[u]’ } {‘-1.6985e-10’} {‘Disease - Reported:Increas…’} {‘Disease:EX_adpoh[u]’ } {‘3.7368’ } {‘Disease - Reported:Increas…’}

2018 - Ines Thiele 2019 - Ines Thiele, included more options (minBiomarker,fixIEMlb,LPSolver) 2021 - Ines Thiele, Tests only for biomarker if the reaction exists in the WBM model

compareBounds2Models(model1, model2)

This function compares the bounds between two models

compareMaleFemale(male, female)

This function compares basic features of the male and female whole-body metabolic models

[ResultsMaleFemale] = compareMaleFemale(male,female)

INPUT male model structure (male whole-body metabolic model) female model structure (female whole-body metabolic model)

OUTPUT ResultsMaleFemale structure containing the basic differences and

commenalities between male and female model

Ines Thiele 2017

convertATPflux2StepNumer(ATP_hydrolysis_flux, sex, weight, height)

This function converts an ATP hydrolysis flux (e.g., Muscle_DM_atp_c_ into distance walked and step number). See below for assumptions and calculation details

function [Energy_kJ,Energy_kcal,Meter,StepNumber] = convertATPflux2StepNumber(ATP_hydrolysis_flux, sex, weight, height) INPUT ATP_hydrolysis_flux Flux value through the Muscle_DM_atp_c_ reaction sex ‘male’ or ‘female’ weight in kg height in cm

OUTPUT Energy_kJ Energy value in kJ corresponding to the Flux value through the Muscle_DM_atp_c_ reaction Energy_kcal Energy value in kcal corresponding to the Flux value through the Muscle_DM_atp_c_ reaction Meter Corresponding meter of walking that can be achieved StepNumber Corresponding step number that can be achieved

Ines Thiele 01/2018

createModelNewCompartment(model, OldComp, NewComp, NewCompName, LB, UB, RemoveExch)

This function converts a two compartment metabolic model into a three compartment metabolic model model

function [modelComp] = createModelNewCompartment(model,OldComp,NewComp,NewCompName,LB,UB,RemoveExch)

INPUT model Model structure OldComp Name of the current (to be replicated) compartment to be

replicated and added

NewCompName Name of new compartment e.g., ‘lumen’ LB Value for Lower bound on new exchange compartment (default:-1000) UB Value for Upper bound on new exchange compartment (default:1000)

NewComp Abrreviation for the new compartment,m e.g., ‘lu’ RemoveExch Option (if set to 1) to remove the old compartment

OUTPUT modelComp Model with three compartments

Ines Thiele 2016-2017

determineFluxValuesOnBoundary(model, solution)

This function determines the number of reactions in the flux distributions that are on the boundaries

[OfConstraint, OfAll] = determineFluxValuesOnBoundary(model, solution)

INPUT model Model structure solution Solution structure

OUTPUT OfConstraint Fraction of flux values that are on the lower and upper bounds

of all constrained reactions (assuming a minimum infinity of -1,000,000 and a maximum infinity of 1,000,000

OfAll Fraction of flux values that are on the lower and upper bounds

of all reactions in the model

Ines Thiele, December 2018

findMetinOrgan(WBModel, metabolite)

find all organs in which a metabolite participates

INPUT WBmodel whole body metabolic model metabolite abbrevition of the metabolite to be looked up

OUTPUT OrganList List of organs that the metabolite occurs in

Ines Thiele, July 2020

getBasicHarveyStats(male, female)

This function compiles basic statistics on the male and female whole-body metabolic model

function [TableHHStats] = getBasicHarveyStats(male,female)

INPUT male Model structure containing the male model female Model structure containing the female model

OUTPUT TableHHStats Table containing the statistics

Ines Thiele 2016- 2019

getListOfUniqueIsozymes(model)

this function gets a (preliminary list of unique isozymes for a given model or reconstruction

INPUT model model structure

OUTPUT listIsozymes list of isozymes

getOrgansFromHarvey(modelWBM, runTestsOnly, OrganCompendium, printLevel)

This function cuts the organs from the whole-body metabolic model. Note that the different biofluid compartments are retained but all constraints on the exchange and transport reactions are overwritten Once the organs are extracted, the function runs the sanity check on each organ. This step can be also done independently by setting runTestsOnly to 1. In this case, the OrganCompendium must be provided as input. The function also loads Recon 3* so that the test results can be compared with the organ test results.

[OrganCompendium,TableCSources] = getOrgansFromHarvey(modelWBM, runTestsOnly, OrganCompendium)

INPUT modelWBM model structure of whole-body metabolic model runTestsOnly

OUTPUT OrganCompendium Structure containing the individal organs as

well as the basic tests that this organ passed

TableCSources Overview table of ATP yield per carbon source

under aerobic and anaerobic conditions for each organ in the model structrure

The organ compendium for each sex will be saved as OrganAtlas_Harvetta.mat and OrganAtlas_Harvey.mat along with the test results.

Ines Thiele, 2017 - 2019

getRxnsFromGene(model, gene, causal)

This function gets all reaction(s) associated with a particular gene by screening through the grRules provided in the model structure

[Rxns, grRules] = getRxnsFromGene(model,gene,causal)

INPUT model model structure gene gene of interest causal if causal == 1 get only genes that would lead to loss of function of the

associated reactions, otw get all associated reactions (default)

OUTPUT Rxns List of reaction(s) associated with the input gene grRules List of grRules containing the input gene, same order as Rxns

Ines Thiele 10/2019

getStatsOrganComp(female, male, OrganCompendium_female, OrganCompendium_male, violinPlots)

This function compiles general statistics on the male and the female organ compendia derived from the male and female whole-body metabolic models.

[TableProp_female,TableProp_male, TableGRM] = getStatsOrganComp(female, male, OrganCompendium_female, OrganCompendium_male, violinPlots)

INPUT female model structure, female whole-body metabolic

model

male model structure, male whole-body metabolic

model

OrganCompendium_female strucutre containing the different organs

(generated with the function getOrgansFromHarvey

OrganCompendium_male strucutre containing the different organs

(generated with the function getOrgansFromHarvey

violinPlots plot violin plots (does not work below Matlab

2016. (defaul = 0)

OUTPUT TableProp_female Table containing organ-specific information TableProp_male Table containing organ-specific information

A more comprehensive comparison output is provided in the file: Results_StatsOrganComp.mat that is created at the end of this function.

Ines Thiele, 2017

linearRegression(x, y)

This function calculates the linear regression in the form of y = a0 + a1*x

function [a0,a1,Rsqr,Residuals] = linearRegression(x,y)

INPUT x Values for the explanatory variable x in y = a0 + a1*x y Values for the dependent variable y in y = a0 + a1*x

OUTPUT a0 Intercept a1 Slope of the line Rsqr Square of the correlation coefficient Residuals Regression residuals, provides an objective measure of the

goodness of fit of the linear regression equation

Ines Thiele 2018

optimizeWBModel(model, param)

Solves flux balance analysis problems, and variants thereof

Solves LP problems of the form

\[\begin{split}max/min ~& c^T v \\ s.t. ~& S v = dxdt ~~~~~~~~~~~:y \\ ~& C v \leq d~~~~~~~~:y \\ ~& lb \leq v \leq ub~~~~:w\end{split}\]

USAGE

solution = optimizeCbModel (model, osenseStr, minNorm, allowLoops, zeroNormApprox)

INPUT

model – (the following fields are required - others can be supplied)

• S - m x 1 Stoichiometric matrix

• c - n x 1 Linear objective coefficients

• lb - n x 1 Lower bounds

• ub - n x 1 Upper bounds

• dxdt - m x 1 change in concentration with time

• csense - m x 1 character array with entries in {L,E,G} (The code is backward compatible with an m + k x 1 csense vector, where k is the number of coupling constraints)

• C - k x n Left hand side of C*v <= d

• d - k x n Right hand side of C*v <= d

• dsense - k x 1 character array with entries in {L,E,G}

OPTIONAL INPUT

• param – Additional parameters as a parameter struct All fields of the struct which are not COBRA parameters (see getCobraSolverParamsOptionsForType) for this problem type will be passed on to the solver in a solver specific manner.

  Some optional parameters which can be passed to the function as part of the options struct (DONE), or as parameter value pairs (TODO), or are listed below:

• * osenseStr – Maximize (‘max’)/minimize (‘min’) (opt, default = ‘max’) linear part of the objective. Nonlinear parts of the objective are always assumed to be minimised.

• * solverName – Solver name {‘tomlab_cplex’,’ibm_cplex’,’cplex_direct’}

• * printLevel – verbose level * if 0, warnings and errors are silenced. (default: 0) * if > 0, warnings and errors are thrown.

• * minNorm – {(0), scalar , n x 1 vector}, where [m, n] = size(S); If not zero then, minimise the Euclidean length of the solution to the LP problem. minNorm ~1e-6 should be high enough for regularisation yet maintain the same value for the linear part of the objective. However, this should be checked on a case by case basis, by optimization with and without regularisation.

OUTPUT

solution – solution object:

• f - Objective value

• v - Reaction rates (Optimal primal variable, legacy FBAsolution.x)

• y - Dual for the molecular species

• w - Reduced costs of the reactions

• s - Slacks of the molecular species

• stat - Solver status in standardized form: * -1 - No solution reported (timelimit, numerical problem etc) * 0 - Infeasible * 1 - Optimal solution * 2 - Unbounded solution

• origStat - Original status returned by the specific solver

• ctrs_y - the duals for the constraints from C

• ctrs_slack - Slacks of the additional constraints

organEssentiality(model, LPSolver)

This function computes the organ essentiality in a whole-body model by setting each organ’s value to zero in the whole-body objective reaction, setting all organ-specific reaction bounds to zero (lower and upper), and then computes whether a non-zero flux through this objective is still possible.

function [ResultsOrganEss] = organEssentiality(model, LPSolver)

INPUT model model structure (whole-body metabolic model) LPSolver Define LPSolver (‘ILOGcomplex’;

‘tomlab_cplex’ -default)

OUTPUT ResultsOrganEss Contains the maximally possible flux value for the

whole-body reaction for each organ. Col 1: Organ name, Col 2: Max flux value, Col 3: Min flux value (if minimization is activated; default: inactive), Col 4: Solver status: 1 = feasible, 5 = feasible with numerical difficulties (rescaling issues), 3 = infeasible)

Ines Thiele 2016 Ines Thiele, added option to specify LPSolver - 10/2018

performIEMAnalysis(model, geneMarkerList, compartment, urine, minRxnsFluxHealthy, causal, reverseDirObj, fractionKO, minBiomarker, fixIEMlb, LPSolver)

This function performs the IEMAnalysis from a list of genes, testing for the defined biomarker metabolites in one or more biofluid compartments.

INPUT model WBM model structure geneMarkerList Cell array containing the geneMarkerLists and

the biomarkers to be tested for e.g., geneMarkerListMarkerList = {

‘5053.1’ ‘trp_L;actyr;phe_L;tyr_L’ ‘249.1’ ‘3pg;cholp;glyc3p;ethamp’ };

compartment List of biofluid compartments that the biomarkers

appear in and should be tested for

urine Indicate whether you want to test for the urine

excretion of the biomarker metabolite as well. Default = true

minRxnsFluxHealthy Min flux value(s) through the IEMRxns (Default: 0.75

corresponding to 75%)

reverseDirObj The function maximizes the objective flux by

default. If set to 1, the function also checks the minimization problem.

fractionKO By default, a complete knowckout of BiomarkerRxnsthe IEM

reactions is computed but it is possible to set a fraction (default = 1 for 100% knockout)

minBiomarker Minimization through biomarker reaction (default = 0) fixIEMlb fix IEM to lb = ub

=(1-fractionKO)*solution.v(find(model.c)) (default = 0, i.e., lb =0, while ub = (1-fractionKO)*solution.v(find(model.c))

LPSolver Define LPSolver (‘ILOGcomplex’ - default;

‘tomlab_cplex’)

causal if causal == 1 get only genes that would lead to loss of function of the

associated reactions, otw get all associated reactions (default)

OUTPUT IEMSolutions Structure containing the predictions for each gene.

Metabolites that are not occurring in a biofluid, will have a ‘NA’ in the corresponding fields

IEMTable Cell array containing the predictions for each gene (same

content as in IEMSolutions

missingMetAll Metabolites not appearing in a biolfuid

Ines Thiele - 2020-2021

performSanityChecksonRecon(model, resultsFileName, ExtraCellCompIn, ExtraCellCompOut, runSingleGeneDeletion, resultsPath, param)

This function performs various quality control and quality assurance tests. [TableChecks, Table_csources, CSourcesTestedRxns, TestSolutionNameOpenSinks,TestSolutionNameClosedSinks] = performSanityChecksonRecon(model,resultsFileName,ExtraCellCompIn,ExtraCellCompOut,runSingleGeneDeletion)

INPUT model model structure resultsFileName File name of the generated output file ExtraCellCompIn [e] compartment by default, if extracellular

uptake compartment is named differently, it can be specified here

ExtraCellCompOut [e] compartment by default, if extracellular

secretion compartment is named differently, it can be specified here

runSingleGeneDeletion if 0 (default): function does not run single gene deletion otw choose 1 resultsPath path where the results should be saved

(default: location where ‘MethodSection3.mlx’ folder is)

param testing parameters. Some test are done by

default, others need to be specified (see below)

OUTPUT TableChecks Table overview of the performed tests and

their outcomes

Table_csources Table the test results for ATP yield from

various carbon scources under aerobic and anaerobic conditions

CSourcesTestedRxns List of reactions active when testing the ATP

yield from the various carbon sources

runIEM_HH

This script predicts known biomarker metabolites in different biofluid compartments (urine, blood, csf) of the whole-body model for 57 inborn-errors of metabolism (IEMs). The meaning of the abbreviations for metabolites and IEMs used in this script can be found at www.vmh.life. The supported model options are ‘male’, ‘female’, and ‘Recon3D’. Please define those using the variable ‘sex’ (e.g., sex = ‘male’).

Ines Thiele 2018 - 2019