Modeling of Metabolic Pathways with MathModelica and Mathematica

Jan Brugård

MathCore Engineering AB

Àlex Gómez-Garrido

Computational Biochemistry and Biophysics Laboratory
Research Unit on Biomedical Informatics
Universitat Pompeu Fabra (UPF)

Preface

To run this notebook you need Mathematica, MathModelica System Designer Professional and the BioChem library.

Please contact your local Wolfram sales representative or MathCore directly to get information about MathModelica System Designer Professional or to get a trial.

The BioChem package will be made available for free from the MathCore web.

To purchase MathModelica System Designer Professional please contact your local Wolfram sales representative or MathCore directly.

Contact Information

MathCore Engineering AB
Teknikringen 1F
583 30 LINKÖPING
Sweden

email: sales@mathcore.com

phone: +46 13 328500

web: www.mathcore.com

The BioBridge Project

1) structured databases including genomic, proteomic, and metabolomic information for the metabolic pathways affected by COPD will be 
	2) integrated to model the underlying metabolic networks. This, in turn, will be utilized by 
	3) a simulation environment for simultaneous analysis of multilevel data able to improve existing knowledge on complex disorders.

Partners:

University of Birmingham
	Institut Für Medizinische und Chemische Labordiagnostik (LKH)
	Genfit
	Biomax Informatics AG
	MathCore Engineering AB
	IDIBAPS
	UPF

Integrative Genomics and Chronic Disease Phenotypes:
modelling and simulation tools for clinicians

Using MathModelica and Mathematica for model development

Development of accurate models require

Estimation of model parameters from measurement data (global optimization, sensitivity analysis)

Based on these experimental data, the parameters of the system can be estimated as they were unknown.

Optimal experimental design can be used to perform further informative experiments and thus to increase the information about the system.

The BioChem Library - A library for modeling of metabolic pathways

General purpose Modelica library for modeling, simulation and visualization of biological and biochemical systems.

The classes implemented in the BioChem library describe substances and reactions that can take place in-between these substances in a diverse number of biochemical pathways.

BioChem provides several packages holding both ready-made models as well as base classes for creating new tailor-made models.

The image below is a screenshow of a pathway developed by, Brands and van Boekel [2]. The model has been imported to MathModelica from the BioModels Database [5].

Building Models in MathModelica

Kinetic model of the Maillard reaction occurring in heated monosaccharide-casein systems imported from BioModels (SBML)

Needs["MathModelica`"];

Simulating

sim1=Simulate[Brands2002⎵MonosaccharideCasein,{t,0,1500}];
PlotSimulation[{Acetic⎵acid.c[t],Amadori.c[t],AMP.c[t],C5.c[t], Fru.c[t],  Glu.c[t]},{t,0,1500}]

Pseudo Experimental Data

Data measured for some of the variables (7 of 11) with 1 minute samples

Load Experiment Data

experimentData=Import["C:\\Users\\brugard\\Desktop\\USA\\Presentation\\Brands2002_60.xls"][[1]];
variableNames={AMP.c[t],Amadori.c[t],Cn.c[t],Formic⎵acid.c[t],Fru.c[t], Glu.c[t],Triose.c[t],lys⎵R.c[t]};
experimentData//MatrixForm

timeVector=experimentData[[1]];
variableValues=Drop[experimentData,1];

ListPlot[variableValues,PlotJoined->True,PlotRange->All]

Parameter Estimation

Using sum of squares

(∑_(i=1)^n)(∑_(j=1)^m)((((Subscript[C, i])[(Subscript[t, j])]-(Subscript[(SuperStar[C]), i])[(Subscript[t, j])]))^2)

where (Subscript[C, i])[(Subscript[t, j])] is the simulated and (Subscript[(SuperStar[C]), i])[(Subscript[t, j])] is the measured concentration of substance i at time (Subscript[t, j])

Residual, alternative 1, Simulate

size=Dimensions[timeVector][[1]];
fResiduals[p1_?NumberQ,p2_?NumberQ,p3_?NumberQ]:=
Module[{sol},
sol=Simulate[Brands2002⎵MonosaccharideCasein,{t,0,1500},ParameterValues->{⎵J1.K1==p1,⎵J2.K2==p2,⎵J3.K3==p3}];
Total[(((Transpose[Table[variableNames/.t->timeVector[[i]],{i,1,size}]]-variableValues))^2),2]
]

fRsim=fResiduals[0.009,0.005,0.0005]//Timing

Residual, alternative 2, NDSolve

ClearAll[$SimulationResult]
{eqns,ieqns,states,vars,pars}=ModelEquations[Brands2002⎵MonosaccharideCasein];
states//MatrixForm

eqns//MatrixForm

initialEqns={
(Melanoidin.n)[0]==0,(AMP.n)[0]==0,(Acetic⎵acid.n)[0]==0,(Cn.n)[0]==0,
(Triose.n)[0]==0,(lys⎵R.n)[0]==15,(Amadori.n)[0]==0,(Glu.n)[0]==160,
(Fru.n)[0]==0,(C5.n)[0]==0,(Formic⎵acid.n)[0]==0,(compartment.V)[0]==1,
(compartment.T)[0]==25
};
variableNames={
(AMP.n)[t],(Amadori.n)[t],(Cn.n)[t],(Formic⎵acid.n)[t],(Fru.n)[t],
(Glu.n)[t],(Triose.n)[t],(lys⎵R.n)[t]
};
fResiduals[p1_?NumberQ,p2_?NumberQ,p3_?NumberQ]:=
Module[{sol},
sol=First[NDSolve[Chop[Join[eqns,initialEqns]/.{⎵J1.K1->p1,⎵J2.K2->p2,⎵J3.K3->p3}/.pars],Join[states,vars],{t,0,1500}]];
Total[(((Transpose[Table[variableNames/.sol/.t->timeVector[[i]],{i,1,size}]]-variableValues))^2),2]
]

Same result as with simulate (but much faster):

fResiduals[0.009,0.005,0.0005]//Timing

fRsim

Finding global minima

Varying k1 and k2, with other variables fixed

ListPlot3D[Flatten[Table[{k1,k2,fResiduals[k1,k2,0.0005]},{k1,0.005,0.02,0.0025},{k2,0.0,0.015,0.001}],1]]

The real problem

In reality all variables K1 to K11 are unknown, and the range is also larger then what we used earlier on.

Find K1 to K11

ClearAll[$SimulationResult]
fResiduals[p1_?NumberQ,p2_?NumberQ,p3_?NumberQ,p4_?NumberQ,p5_?NumberQ,p6_?NumberQ,p7_?NumberQ,p8_?NumberQ,p9_?NumberQ,p10_?NumberQ,p11_?NumberQ]:=
Module[{sol},
sol=First[NDSolve[Chop[Join[eqns,initialEqns]/.{⎵J1.K1->p1,⎵J2.K2->p2,⎵J3.K3->p3,⎵J4.K1->p4,⎵J5.K5->p5,⎵J6.K6->p6,⎵J7.K7->p7,⎵J8.K8->p8,⎵J9.K9->p9,⎵J10.K10->p10,⎵J11.K11->p11}/.pars],Join[states,vars],{t,0,1500}]];
Total[(((Transpose[Table[variableNames/.sol/.t->timeVector[[i]],{i,1,size}]]-variableValues))^2),2]
]

NMinimize[
Max[0,fResiduals[k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11]],
{{k1,0.001,0.1},
{k2,0.0001,0.001},
{k3,0.00001,0.001},
{k4,0.0001,0.001},
{k5,0.0001,0.01},
{k6,0.0001,0.01},
{k7,0.00001,0.001},
{k8,0.01,1},
{k9,0.01,0.1},
{k10,0.00001,0.001},
{k11,0.01,1}
},
AccuracyGoal->4,PrecisionGoal->3,MaxIterations->50000,
StepMonitor:>Print["Parameters: K1 = ",k1,"     K2 = ",k2,"     K3 = ",k3]]

solution=%[[2]]

Simulating Again

sim2=Simulate[Brands2002⎵MonosaccharideCasein,{t,0,1500},ParameterValues->{⎵J1.K1==k1,⎵J2.K2==k2,⎵J3.K3==k3,⎵J4.K4==k4,⎵J5.K5==k5,⎵J6.K6==k6,⎵J7.K7==k7,⎵J8.K8==k8,⎵J9.K9==k9,⎵J10.K10==k10,⎵J11.K11==k11}/.solution];

PlotSimulation[{Acetic⎵acid.c[t], Fru.c[t],  Glu.c[t]},{t,0,1500},SimulationResult->{sim1,sim2}]

Optimal Experimental Design - Design Criterias

F = Fisher information matrix
∑ = (F^(-1)) = Covariance matrix

A-optimal design, max(Tr(F))

This criterion is rarely used since it can lead to non informative experiments with a covariance matrix that is not positive definite.

D-optimal design, max(det(∑))

It can be interpreted as minimizing the geometric mean of the errors in the parameters. It consists in minimizing the total parameter uncertainty in the
system.

E-optimal design, min((Subscript[λ, max])(∑	))

It minimizes the largest error, which corresponds to a minimization of the largest eigenvalue of the covariance matrix. It consists in minimizing the
parameter uncertainty in the least constrained direction in the parameter space.

Modified E-optimal design, min((((Subscript[λ, max])(∑))/((Subscript[λ, min])(∑))))

As it minimizes the ratio of the largest to the smallest eigenvalue, it optimizes the functional shape of the confidence intervals.

Optimal Experimental Design, Covariance Matrix

Parameters to fit

kineticParameters = {⎵J1.K1, ⎵J2.K2, ⎵J3.K3, ⎵J4.K4, ⎵J5.K5, ⎵J6.K6, ⎵J7.K7, ⎵J8.K8, ⎵J9.K9, ⎵J10.K10, ⎵J11.K11};
averageConcentrations=(1/1500){NIntegrate[AMP.c[τ],{τ,0,1500}],NIntegrate[Acetic⎵acid.c[τ],{τ,0,1500}],NIntegrate[Amadori.c[τ],{τ,0,1500}],NIntegrate[C5.c[τ],{τ,0,1500}],NIntegrate[Cn.c[τ],{τ,0,1500}],NIntegrate[Formic⎵acid.c[τ],{τ,0,1500}],NIntegrate[Fru.c[τ],{τ,0,1500}], NIntegrate[Glu.c[τ],{τ,0,1500}],NIntegrate[Melanoidin.c[τ],{τ,0,1500}],NIntegrate[Triose.c[τ],{τ,0,1500}],NIntegrate[lys⎵R.c[τ],{τ,0,1500}]};
nConc[t_]:={AMP.c[t],Acetic⎵acid.c[t],Amadori.c[t],C5.c[t],Cn.c[t],Formic⎵acid.c[t],Fru.c[t], Glu.c[t],Melanoidin.c[t],Triose.c[t],lys⎵R.c[t]}/averageConcentrations;

Calculate sensitivity

Simulate[Brands2002⎵MonosaccharideCasein , {t, 0, 1500}];
originalConcentrations=Table[nConc[t],{t,0,1500,1.5}];
sensitivityCalculator[i_]:=Module[
{deltaConcentrations,ΔT=1.5,δ=(10^(-6))},
Simulate[Brands2002⎵MonosaccharideCasein , {t, 0, 1500}, ParameterValues -> { kineticParameters[[i]]==  (1+δ)(kineticParameters[[i]] /. pars)}];
deltaConcentrations = Table[nConc[t],{t,0,1500,ΔT}]; 
((Abs[(deltaConcentrations - originalConcentrations)])/(δ kineticParameters[[i]] /. pars)) 
];

sensitivityTable=Table[sensitivityCalculator[i],{i,1,11}];

Fisher Matrix and Covariance

getFIMcoefficient[a_, b_]:= Total[sensitivityTable[[a]]sensitivityTable[[b]],2]
FIM = Table[getFIMcoefficient[a,b], {a, 1, 11}, {b, 1, 11}] ;
COV = Inverse[FIM];
CorrelationMatrix = Table[((COV[[i]][[j]])/(Sqrt[COV[[i]][[i]] COV[[j]][[j]]])), {i, 1, 11}, {j, 1, 11}] ;

MatrixPlot[CorrelationMatrix]

Summary

New Library

A new library, BioChem, for modeling and simulation of metabolic pathways will soon be available for free

SBML to Modelica converter

An SBML to Modelica converter is being developed. This will allow users to import existing SBML models to Modelica

Analysis

The integration with Mathematica allows users to perform both parameter estimation and optimal experimental design

Results from MathModelica can be converted into Mathematica equations and used by Mathematica functions such as NDSolve (If ODE).

COPD

The results of the BioBridge project will be used to further understand COPD.

Acknowledgements and References

Acknowledgements

The first version of the library was created by Emma Larsdotter Nilsson and Erik Ulfhielm at Linköping University. The current version of the library has been further developed by  MathCore.

The BioBridge project is supported by the European Commission and is part of the Sixth Framework Programme

References

1. Leonardo M Fabbri, Klaus F Rabe, From COPD to chronic systemic inflammatory syndrome?, Lancet 2007; 370: 797–99
2. Brands CM , van Boekel MA, Kinetic modeling of reactions in heated monosaccharide-casein systems, J. Agric. Food Chem. (ISSN: 0021-8561), November 2002
3. Gomez-Garrido, A. López García de Lomana, J. Villà-Freixa, Parameter Estimation for Biochemical Models in a Variety of Parameter Landscapes, Bioinformatics 2007. In preparation 
4. A. López García de Lomana, A Gomez-Garrido, M. Hernández-Sanchez, J. Villà-Freixa, ByoDyn: an integrative tool for computational systems biology, Bioinformatics 2007. Submitted

5. BioModels Database, http://www.ebi.ac.uk/biomodels/
6. BioBridge, http://www.biobridge.eu/ 
7. MathCore Engineering AB, http://www.mathcore.com/, info@mathcore.com