Analysis of complex reaction networks using mathematical programming - - PowerPoint PPT Presentation

Marianthi Ierapetritou

Department Chemical and Biochemical Engineering Piscataway, NJ 08854-8058

Analysis of complex reaction networks using mathematical programming approaches

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Complex Process Engineering Systems?

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

General Motivation

Diverse complex systems spanning different scales

Liver metabolism (molecular level) Combustion systems (process level) Scheduling of multiproduct-multipurpose plants (plant

level)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Motivation -1: Liver Support Devices

Acute and chronic liver failure account for 30,000

deaths each year in the US

A large number of liver diseases:

• Alagille Syndrome
• Alpha 1 - Antitrypsin Deficiency
• Autoimmune Hepatitis
• Biliary Atresia
• Chronic Hepatitis
• Cancer of the Liver
• Cirrhosis
• Cystic Disease of the Liver
• Fatty Liver
• Galactosemia
• Hepatitis A, B, C

Currently liver transplantation is primary therapeutic

• ption. Scarcity of donor organs limits this treatment

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Solutions

Adjunct Internal Liver Support With Implantable

Devices

Hepatocyte Transplantation Implantable Devices Encapsulated Hepatocytes

Extracorporeal Temporary Liver Support

Nonbiological devices: hemodialysis, hemofiltration, plasma exchange units Hepatocyte- and liver cell–based extracorporeal devices

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

a) How to maximize long-term functional stability of hepatocytes in inhospitable environments b) How to manufacture a liver functional unit that is scalable without creating transport limitations or excessive priming volume that must be filled by blood

• r plasma from the patient

c) How to procure the large number of cells that is needed for a clinically effective device

Challenges

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

… and the Reality

Problem complexity: System of large interconnectivity Large number of adjustable variables Uncertainty

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Accuracy depends on :

• Flow model
• Kinetic model

Fluid flows significantly affected by chemical reaction : Combustion, Aerospace propulsion

Conversion of chemical energy to mechanical energy Require alternate representation of complex kinetic mechanism, without sacrificing accuracy

Motivation – 2 : Combustion

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Challenges: Combine Flow and Chemistry

How should these be combined ?

Composition map Velocity map

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

…and the Reality

Complex kinetics

(LLNL Report, 2000)

H2 mole fraction vs. time

Uncertainty in kinetic parameters

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Motivation-3: Large-Scale Process Operations

Crude Oil Marine Vessels Storage Tanks Charging Tanks Crude Distillation Units Other Production Units Component Stock Tanks Blend Header Finished Product Tanks Lifting/ Shipping Points

Product Blending & Distribution Crude-oil Unloading and Blending Production

Max Profit

Subject to: Material Balance Constraints Allocation Constraints, Sequence Constraints Duration Constraints, Demand Constraints …

Goal: Address the optimization of large-scale short-term scheduling problem, specifically in the area of refinery operations

Add $100s of million/year profit by

• ptimizing crude-oil-

marketing enterprise

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Challenges: Parameter Fluctuations

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

10% 90% Product 2 Product 1 Feed A Hot A 40% 60% 60% Int AB 80% Feed B 50% Feed C Imp E 50% 20% Int BC Heating Reaction 2 Reaction 1 Reaction 3 Separatio n 40%

T Ti im me e N Nu um mb be er r

• f

f E Ev ve en nt t P Po

• i

in nt ts s O Ob bj je ec ct ti iv ve e f fu un nc ct ti io

• n

n v va al lu ue e C CP PU U t ti im me e ( (s se ec c) ) 8 8 h ho

• u

ur rs s 5 5 1 14 49 98 8. .1 19 9 0. .4 47 7 1 16 6h ho

• u

ur rs s 9 9 3 37 73 37 7. .1 10 1 17 77 7. .9 93 3 2 24 4 h ho

• u

ur rs s 1 13 3 6 60 03 34 4. .9 92 2 9 92 23 36 67 7. .9 94 4

As time horizon of scheduling problem increases, the solution requires exponential computational time which makes the problem intractable.

…and the Reality

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Systems Approaches

Mathematical programming

Systematic consideration of variable dependences Continuous and discrete representation

Sensitivity – parametric analysis

Identification of important features and parameters

Feasibility evaluation

Conditions of acceptable operation

Optimization

Multiobjective since we have more than one objective to

• ptimize

Uncertainty

Evaluation of solutions that are robust to highly fluctuating environment

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Presentation Outline

Complexity reduction using mathematical

programming approaches

Optimization of hepatocyte functionality Reduction of complex chemistry Uncertainty analysis & feasibility evaluation Analysis of alternative solutions

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Hepatic Metabolic Network

45 internal metabolites 76 reactions: 33 irreversible + 43 reversible 34 measured (red) + 42 unknown

Chan et al (2003) Biotechnol & Bioengineering

Main Assumptions

1) Gluconeogenic and fatty acid oxidation enzymes are active in plasma 2) Energy-requiring pathways are negligible 3) Metabolic pools are at pseudo-steady state.

Main Reactions

Glucose Metabolism (v1-v7) Lactate Metabolites & TCA Cycle(v8-v14 ) Urea Cycle (v15-v20) Amino acid uptake & metabolism (v21-v68, ,v76) Lipid & Fatty Acid Metabolism (v46-v50,v71-v75)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Measure 2 fluxes: Uniquely-determined system Measure 3 fluxes: Overdetermined System- Least Square method Measure 1 flux: Underdetermined System- Linear Programming

Pseudo-steady State

= 0 Sv

1 2 3 4 5 1 2 3 4

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 dA v dt v dB dt v dC dt v dD v dt b dE dt b dF b dt dN b dt ⎡ ⎤ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ − − ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ = − − ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

A B C D E F N 2 N v1 v2 v3 v4 v5 b1 b2 b3 b4

= −

u u m m

S v S v

Metabolic Flux Analysis (MFA) is developed to calculate unknown intracellular fluxes based on the extracellular measured fluxes.

Rationale for Metabolic Modeling

Interpretation and coupling to experimental data. Gain insights into how cells adapt to environmental changes. To identify key pathways for hepatocyte function.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Optimization in Metabolic Networks

• Single-level Optimization: Optimize a single objective function

(e.g. maximization of a single metabolic flux).

Multi-Objective Optimization Multi-level Optimization

≠

Uygun et al., (2006) Ind. Eng. Chem. Res. Lee S. et al (2000) Computer & Chem. Eng.

• Schilling. et al., (2001) Biotechnol Bioeng

Eward & Palsson (2000) PNAS Segre D. et al (2002) PNAS

• Multi-level Optimization: Several objectives acting hierarchically to optimize

their own objective function (e.g. Minimize the difference of predicted fluxes from experimentally observed values to optimize the cellular objective function).

Nolan R.P et al (2005) Metabolic Engineering Burgard & Maranas (2003) Biotechnol Bioeng Uygun et al., (2007) Biotechnol Bioeng Sharma N.P. et al., (2005) Biotechnol Bioeng

• Multi-objective Optimization: Several objective functions are

simultaneously optimized (e.g. minimizing the toxicity and maximizing metabolic production).

Nagrath D. et al. (2007) Annals of Biomedical Engineering

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Aim: Identify the flux distributions for optimal urea production that can fulfill metabolites balances and flux constraints

Unit: µmol/million cells/day

1 min max

: :

urea N ij j j j j j

Max Z v Subject to S v i M v v v j K

=

= = ∀ ∈ ≤ ≤ ∀ ∈

∑

Single-level Optimization: Maximize Urea Secretion

> 2 fold 2.35±0.52 LPAA > 15 fold 0.17±0.24 LIP > 3 fold 1.32±0.69 HPAA > 10 fold 6.81 0.23±0.43 HIP Increase Optimal Value Experimental Data*

*Chan & Yarmush et al (2003) Biotechnol Prog

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Results for Single-Level Optimization

Increased fluxes Gluconeogenesis (R2-R6) TCA Cycle (R13,R14) Urea Cycle (R16,R17) Amino Acid Catabolism (R21,R23,R27,R30,R36,R38,R43) Fatty Acid Metabolism (R47,R48) Pentose Phosphate Pathway (R54) Amino acid uptake fluxes (e.g: Arginine, Serine, Glycine.....)

Fluxes significantly altered through the pathways (more than 30 % change)

Decreased fluxes Amino Acid Catabolism (R33,34) Fatty Acid Oxidation (R46) Glycerol uptake and metabolism, glycogen storage (R70,R71,R73,R74)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Multi-objective Optimization

Sharma NS, Ierapetritou MG,Yarmush ML.,Biotechnol Bioeng. 2005 Nov 5;92(3): 321-35.

albumin N ij j i j 1 min max j j j urea

max v subject to : S v b, i 1,...., M v v v v ε

=

= = ≤ ≤ ≥

∑

• Pareto optimal set yields the feasible

region for BAL operation

• Point A,D and B belong to the Pareto Set
• Both urea and albumin secretion can be

improved at point C by moving towards the Pareto set

Results ε-Constraint Method

Different values of ε were used to calculate the maximum albumin Production = Pareto set

0.0001 0.0002 0.0003 0.0004 0.0005 6.69577 6.71836 6.74095 6.76353 6.78612 6.80871 Urea Secrection (µmol/million cells/day) Albumin Secrection (µmol/million cells/day)

C D A B Feasible Region

0.0001 0.0002 0.0003 0.0004 0.0005 6.69577 6.71836 6.74095 6.76353 6.78612 6.80871 Urea Secrection (µmol/million cells/day) Albumin Secrection (µmol/million cells/day)

C D A B Feasible Region

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina Modulation by optimal AA. Supplementation

Optimal Amino Acid Supplementation

0.5 1 1.5

P r

• l

i n e S e r i n e H i s t i d i n e G l y c i n e A l a n i n e A r g i n i n e T y r

• s

i n e I s

• l

e u c i n e L e u c i n e T h r e

• n

i n e V a l i n e

Concentration (m M)

Optimal Flux Distribution

Culture media supplementation to improve

cellular function

Advantage is no direct genetic intervention

Assumption: Linear Relationship

Higher Amino Acid Supplementation Similar & Low Amino Acid Supplementation Plasma AA supplementation in low-insulin Optimal AA supplementation

1 2 3 4 5 6 7

A p a r t a t e A s p a r a g i n e C y s t e i n e G l u t a m a t e M e t h i

• n

i n e P h e n y l a l a n i n e L y s i n e G l u t a m i n e

Concentration (mM)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Bi-level Optimization

Leader Objective: Minimize the difference between the native type and the knockout condition after reaction deletion

NA j

v : represents the flux distribution

• f native type determined from MILP

model

Aim: Compare the flux distributions between the wild-type and

knock out condition and identify the essential reactions for target cell functions

Segre D. et al (2002) PNAS; 99: 15112-15117 Burgard et al., (2003) Biotechnol Bioeng; 84: 647-657

2 , 1 , . . : . . | | :

max min 1

= ≤ ≤ = = −

∑ ∑

= ∈ d j j j N j j ij urea N j NA j j

v v v v M i v S t s v Max t s v v Min … … Follower Objective: Maximize the particular cell function (urea production)

Critical Pathways for Urea and Albumin Function

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

R16,17 R43 R2,3,4 R5,6 R36 R38 R30 R47 R54 R73 R7 R39 R23 R32 R48 R21 R34 R50 R25 R57 R8 R40 R55 R33 R56 R45 R74 R71 R69 R72 R1 R75 0.2 0.4 0.6 0.8 1 Deleted Reaction Zknockout/Z

Important Pathways

• ptimal value

after individual reaction deletion

• ptimal value of

native type

Gluconeogenesis(R2-R6), Urea Cycle (R16,R17)

:

knockout

Z

: Z

Amino Acid Catabolism (R30,R36,R38,R43) Fatty Acid Metabolism (R47), Pentose Phosphate Pathway (R54)

Using Extended Kuhn-Tucker Approach

• C. Shi. (2005) Applied Mathematics & Computation

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Comparison of Different Methodologies

II) AA supplements Low Insulin (50 µU/ml) III) no supplements IV) AA supplements I) no supplements High Insulin (0.5 U/ml)

Experiments: Model:

Subject to:

Results:

Subject to:

6.809 59.886 Primal-Dual 2.254 0.246 KKT Urea Error Approach Case 2:

NA j

v

• Measured fluxes from Experiment LPAA

6.809 269.239 Primal-Dual 0.165 0.079 KKT Urea Error Approach Case 1:

NA j

v

• Measured fluxes from Experiment LIP

k j v v v M i v S v Max v v Min

j j j N J j ij urea fluxes Measured j NA j j

∈ ∀ ≤ ≤ = = −

∑ ∑

= ⋅ ∈

, 2 , 1 , : | | :

max min 1

… …

Hong, Roth, Ierapetritou, AIChE Annual Meeting, Nov 2007

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Critical Pathways for Urea and Albumin Function

j

N j j 1 N ij j i j 1 m in m ax j j j j j

m in subject to: S v b , i 1, ..., M v v v , j 1, ..., N

λ = =

Φ = λ = = λ ≤ ≤ λ =

∑ ∑

Logic based programming

λj is a binary variable corresponding to the presence or absence of reaction (j) in the network.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Critical pathways for urea and albumin function

Unsupplemented High Insulin HIP Optimal Amino Acids Low Insulin Optimal Amino Acids Low Insulin LIPAA Unsupplemented Low Insulin LIP Plasma Supplementation Medium Pre-Conditioning Condition

Different Conditions

Elucidate Insulin Effects Elucidate AA Effects

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Optimal Condition vs. LIPAA

• Compensatory effects in

TCA cycle fluxes

• Lower gluconeogenic and

lipid metabolism pathway fluxes.

• Higher urea cycle fluxes.
• Higher AA uptake rates

Thick red lines correspond to higher fluxes for optimal condition as compared to LIPAA. Thick blue lines correspond to lower fluxes for optimal condition as compared to LIPAA. Dotted red lines correspond to reactions not important in Optimal case for maximal urea and albumin function.

Glycogen Glucose-6-P Glucose 2-G3P PEP Pyruvate Lactate Glycerol-3-P PPP Glycerol Triglyceride storage Triglyceride Tyrosine Leucine Tyrosine Leucine Ketone Bodies ACAC Phenylalanine Phenylalanine Arginine Aspartate Ac-CoA OAA Cysteine Serine Malate Cysteine Serine Threonine Alanine Glycine Prop-CoA Tyrosine Valine Methionine Isoleucine Glutamine Histidine Proline Glutamate Citrate Fumarate Alpha- ketoglutarate SCC-CoA 14 Glycine NH3 61,63, 62 Glutamine Histidine Proline Lysine ACAC-CoA Citrulline Ornithine Arginine NH3 Ornithine 20 Urea Fatty Acids Fatty Acids CoA 46 5 2,3,4 Aspartate 43 Asparagine Asparagine NH3 NH3 Respiratory Chain NADH H2O NAD+ O2 51 17 15 18 16 16 19 Valine Methionine Isoleucine 64,66,67 11 12 13 6 7 Methionine 24 Alanine 22 8 34 29 28 Glutamate 37 36 10 9 Lysine Threonine O2 53 39 23 9 21 25 31 32 48 34 76 71 54 72 42 38,41,40 4 8 60 55,42,56 57 4 7 35,68 3 3 59 58 65 45 70 74 73 1 27 Tryptophan 30 ADP +Pi ATP 26 4 4 4 9 , 5

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Complexity reduction using mathematical

programming approaches

Optimization of hepatocyte functionality Reduction of complex chemistry Uncertainty analysis & feasibility evaluation Analysis of alternative solutions

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

H2 + O2 2OH OH + H2 H2O+H O + OH O2+H O + H2 OH+H H + O2 HO2 OH + HO2 H2O+O2 H + HO2 2OH O + HO2 O2+OH 2OH O+H2O H + H H2 H + H + H2 H2 + H2 H + H + H2O H2 + H2O H + OH H2O H + O OH O + O O2 H + HO2 H2 + O2 HO2 + HO2 H2O2 + O2 H2O2 OH + OH H2O2 + H HO2 + H2 H2O2 + OH H2O + HO2

Detailed kinetic models are extremely complex Detailed kinetic models are extremely complex

Reduction of complex kinetic mechanism to enable detailed flame simulation

Model Reduction Using Mathematical Programming

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Objective function :

i N N i

R S

∑

= / 1

λ

min Constraint : δ χ ≤

represents total number of species / reactions

i N N i

R S

∑

= / 1

λ

where χ is an error measure representing deviation of full profile from reduced profile

i

λ

Constraint : retain desired system behavior within prescribed accuracy : Binary variable corresponding to ith reaction/species

Optimization Based Reduction

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

s k k k

N k M R dt t dy , , 1 ) ( … = ρ =

∑

=

− =

s

N k p k k k

C h M R dt dT

1 ρ i f ki r ki N i k

q R

r i

) (

1

γ − γ λ =∑

=

(Solved using DVODE) Detailed mechanism Reduced model Time (sec) Temperature (K) Detailed mechanism Reduced model Time (sec) CH4 mass frac

Perfectly Stirred R.

(absolute mixing)

Perfectly Stirred R.

(absolute mixing)

Constraint : δ χ ≤

CH4=0.055, O2=0.19 T=1200 K Methane mechanism(GRI-3.0) : 53 species, 325 reactions Reduced Model : 17 species, 59 reactions, δ = 0.085

Initial Condition:

2 1

2 2

) ( ) ( ) ( ) ( ) ( ) ( ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − =

∑ ∫ ∫

∈κ

χ

k t t full full reduced full k full k reduced k

dt t T t T t T dt t y t y t y

δ χ ≤ Constraint :

Batch Reactor

Evaluation of Constraint Function

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Perform reaction reduction with

Nr binary variables

Two Step Solution Procedure

Binary variables for species reduction (Ns) : 53 Binary variables for reaction reduction (Nr) : 325 Species reduction Eliminate reactions associated with removed species Generate initial reduced reaction set (Ns < Nr)

Final reduced model Mathematical Model: MINLP with embedded ODEs

Methane mechanism: GRI 3.0 (17 species, 113 reactions) (17 species,59 reactions)

Banerjee and Ierapetritou, Chem. Eng. Sci, 8, 4537, 2003.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Detailed model Reduced model Detailed model Reduced model Watched Species Non-Watched Species Reduced Model : 17 species, 59 reactions, δ = 0.085 CH4=0.26, O2=0.086 T=1200 K Reduced model Detailed model

Temperature (K ) CO mass frac. CH3 mass frac. Time (s)

Detailed model Reduced model

CH4 mass frac. Time (s) Time (s) Time (s)

Performance of Reduced Models

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Ns Nr Sparsity Cummulative(CPU) 53 325 1227 (0.07) 190.3 29 126 461 (0.126) 29.57 22 81 291 (0.163 13.87 22 35 131 (0.17) 5.67 19 59 210 (0.187) 5.75 20 30 112(0.187) 3.9 20 25 95 (0.19) 2.34 20 22 84 (0.19) 1.83 Full Mechanism Reduced Mechanisms

Computational Savings by Reduction

50% 96%

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Detailed mechanism Reduced model

Nominal condition

Temperature (K ) Time (s)

Reduced model Detailed mechanism

Outside feasible range

Temperature (K ) Time (s)

Reduced Model has Limited Range of Validity

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Complexity reduction using mathematical

programming approaches

Optimization of hepatocyte functionality Reduction of complex chemistry Uncertainty analysis & feasibility evaluation Analysis of alternative solutions

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Feasibility Quantification

Pressure Temperature Safe Operating Regime

Determine the range operating conditions for safe and productive operations Given a design/plant or process

Design

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Feasible Range Desired Range of Variability

Feasibility Quantification

Convex Hull Approach

(Ierapetritou, AIChE J., 47, 1407, 2001)

Systematic Way of Boundary Approximation

Flexibility Range (Grossmann and coworkers) Deviation of nominal conditions

Nominal Value of Product 1 Nominal Value

• f Product 2

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

A powerful, approach available to identify the uncertainty ranges where the design, process or material is feasible to operate or function.

Process Flexibility

θ1 θ1

N

θ2

N

θ2

) , ( = θ ψ d

∆Θ1+ ∆Θ1

• ∆Θ2+

∆Θ2-

+ +

Δ − =

j N j j j

θ θ θ δ

− −

Δ − =

j j N j j

θ θ θ δ p j ,.., 1 =

T (Swaney & Grossmann 1985)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

δ1 δ2 F

Flexibility Index

+ −

Δ + ≤ ≤ Δ − θ θ θ θ θ F F

N N

Feasible operation can be guaranteed for T Flexibility Index F – one-half the length of the side of hypercube T

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

δ min = FI

. .t s

) , ( = θ ψ d

u d

z

min ) , ( = θ ψ

I i x z d hi ∈ = , ) , , , ( θ

( , , , ) ,

j

g d z x u j J θ ≤ ∈

} | { ) (

+ −

Δ − ≤ ≤ Δ − = θ δ θ θ θ δ θ θ δ

N N

T

≥ δ

Mathematical Formulation

(Swaney & Grossmann 1985)

( )

max minmax

z T i I θ δ ∈ ∈

) , , ( ≤ θ z d fi

} | { ) (

+ −

Δ − ≤ ≤ Δ − = θ δ θ θ θ δ θ θ δ

N N

T

. .t s

δ max = FI

Feasibility test

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Active Set Strategy

δ min = FI

. .t s

I i x z d hi ,..., 1 , ) , , , ( = = θ

( , , , ) 0, 1,...,

j j

g d z x s j J θ + = =

1

1

=

∑

= J j j

λ

1 1

= ∂ ∂ + ∂ ∂

∑ ∑

= = I i i i J j j j

z h z g μ λ

≤ −

j j

w λ

) 1 ( ≤ − −

j j

w U s

1

1

+ =

∑

= z J j j

n w

+ −

Δ − ≤ ≤ Δ − θ δ θ θ θ δ θ

N N

}, 1 , { =

j

w , , ≥

j j s

λ δ

Inner problem is replaced by KKT constraints (Grossmann & Floudas 1987)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Simplicial Approximation

1 2 3 Find mid point of largest tangent plane Insert the largest hypersphere in the convex hull Find Convex hull with these points (1-2-3) Choose m n+1 points for n dimensions (points 1,2,3) 1 2 3

≥

1 2 3 Find new boundary points by line search from the mid point 4 Inflate the convex hull using all the new points After 4 iterations Approximate Feasible Region1- 2-3-4-5-6-7

Goyal and Ierapetritou, Comput. Chem. Engng. 28, 1771, 2004

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Noncovex Problems: Need for Alternative Methods

Failure of Existing Methods due to Convexity Assumptions Assumption: The Non- Convex Constraints can be identified a priori

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Alpha-shape method: Eliminate maximum possible circles of radius α without eliminating any data point For the α shape degenerates to the

• riginal point set

For the α shape is the convex hull of the original point set

(Ken Clarkson http://bell-labs.com/netlib/voronoi/hull.html)

Improved Feasibility Analysis: Shape Recosntruction

α →

α → ∞

Given a set of points, determine the shape formed by these points Analogous to problem of shape reconstruction Problem definition : Given a set of points (sample feasible points), determine mathematical representation of occupied space

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Disjoint nonconvex object Conventional techniques of inscribing hyper-rectangle or convex hull performs poorly

Inscribed convex hull

Boundary points identified by α shape

Identify boundary points using α shape

Polygonal approximation

Connect boundary points to form a polygon

Improved Feasibility Analysis by α - Shapes

Banerjee and Ierapetritou, Ind. Eng. Chem. Res., 44, 3638, 2005.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Broad range of species concentration and temperature encountered in flow simulation Different reduced models for different conditions encountered in flow simulation

Reduced Set # 1 Reduced Set # 2 Reduced Set # 3 Reduced Set # 4

Adaptive Reduction

Banerjee and Ierapetritou, Comb. Flame, 144, 219, 2006.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Sample the feasible space Construct α – shape with the sampled points Determine points forming the boundary of the feasible region

Estimation of Feasible Region: α –shape

20 reduced sets

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Model 1 Species = 22, Reactions = 81 Information of feasible region Model 2 Species = 19, Reactions = 59 Information of feasible region Model 3 Species = 20, Reactions = 22 Information of feasible region : : Model 20 Species = 53, Reactions = 325 Information of feasible region

Library of reduced models

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + + = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ =

s s s s

n n n t n t

S uy uy uy u P E P u u F y y y E u W ρω λ ρω λ ρω λ ρ ρ ρ ρ ρ ρ ρ ρ ρ ρ

• 2

2 1 1 2 1 2 2 1

, ) ( ,

y1,y2, …, yns, T (Checks for a feasible model) Determine λ1, λ2, …, λns Set 6 Set 7 Set 3 Set 1 Set 5 Set 4 Set 2

, S x F t W = ∂ ∂ + ∂ ∂

Reactive flow model

Generate Library of Reduced Model

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina Detailed model Adaptive Reduced model Single Reduced model Detailed model Adaptive Reduced model Single Reduced model Detailed model Detailed model Single Reduced model Single Reduced model Adaptive Reduced model Adaptive Reduced model

Temperature (K ) Time (s) Time (s) Time (s) Time (s) CH4 mass frac. H mass frac. H2 mass frac.

Single reduced model : 38% error Adaptive reduced model : 3% error

Adaptive Reduction Model in PMSR Simulation

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Uncertainty in kinetic parameters

Uncertainty inherent in kinetic parameter data Commonly characterized by

Error bounds (Δlogkf,i, ΔEi etc.), confidence intervals/ranges. Multiplicative Uncertainty Factor (UF ≥1) Upper bound = UF*kf,i, Lower Bound = kf,i/UF Objective: Development of an accurate, systematic and efficient framework of analysis, that characterizes uncertainty in kinetic mechanisms

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − =

β

T R E T A k

g i a i i f

i

, ,

exp

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Representation of Uncertainty

Classical/Rough Set Theory, Fuzzy Measure/Set Theory,

Interval Mathematics and

Probabilistic/Statistical Analysis

Sensitivity Testing Methods Analytical Methods – Differential Analysis e.g. Perturbation Methods – Green’s Function Method – Spectral Based Stochastic Finite Element Method forms the basis of the Stochastic Response Surface Method (SRSM) Sampling Based Methods e.g. – Monte Carlo Methods – Latin Hypercube Methods

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Stochastic Response Surface Method

• Extension of classical deterministic Response Surface

Method and newer Deterministic Equivalent Modeling Method

• The outputs are represented as a polynomial chaos

expansion in terms of Hermite polynomials:

• Allows for direct and probabilistic evaluation of statistical

parameters of the outputs e.g., for the second order

• utput U2: Mean = α0,2

Variance =

∑

=

ξ + =

n i i i

a a U

1 1 , 1 , 1

∑∑ ∑ ∑

− = > = =

ξ ξ + − ξ + ξ + =

1 1 2 , 1 2 2 , 1 2 , 2 , 2

) 1 (

n i j i ij n i j n i i ii n i i i

a a a a U

2 2 1,2 11,2

2 a a +

1st order 2nd order

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Stochastic Response Surface Method

• Method Outline:
• Choice of order of expansion and transformation of the

set of parametric input uncertainties in terms of a set of standard random variables (srv’s) ξ’s - Gaussian (N(0,1)). Commonly encountered transformations include :

Exponential Lognormal (μ,σ) Normal (μ,σ) Uniform (a,b) Transformation Distribution Type

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ξ + λ − ) 2 / ( 2 1 2 1 log 1 erf

( )

σξ + μ exp

σξ + μ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ξ + − + ) 2 / ( 2 1 2 1 ) ( erf a b a

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Stochastic Response Surface Method

• Generation of input points following the Efficient

Collocation Method (ECM)

• Points are selected from the roots of Hermite

polynomials of higher order than the expansion

• Borrows from Gaussian quadrature
• Application of the model to these input points and

computation of relevant model outputs

• Estimation of the unknown coefficients of the expansion

via regression using singular value decomposition (SVD)

• Statistical and direct analysis of the series expression
• f the outputs

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

SRSM - Algorithm

Input Distributions Select a set of srv’s and transform inputs in terms of these Express outputs as a series (of chosen order) in srv’s with unknown coefficients Generate a set of regression points

Model

Estimate the coefficients of the

• utput approximation

Output Distributions

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Implementation

Discretization of time interval 2nd order SRSM expansion fit for each output

species at each time point

MODEL

y = f(k1, k2, …, k19)

A19 P(A19)

. . .

P(y) y A1 P(A1)

Initial Conditions

}

P(y) y

y t t=t1 t=t2

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Uncertainty Propagation: Results

Concentration profiles display time varying distributions Number of model simulations required by SRSM is orders of

magnitude less than Monte Carlo (723 vs. 15,000)

Distributions of H2 at t=1, 2 and 5 seconds generated by Monte Carlo (MC) simulation and SRSM Nominal H2 mole fraction vs. time plot

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Uncertainty Propagation: Results

Output distributions at

each time point very well approximated by second

• rder SRSM

Sensitivity information

easily obtained via expansion coefficients - aids understanding how the reaction sequence progresses

Means for successfully

preprocessing the reduction of the kinetic model taking into account uncertainty

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Presentation Outline

Complexity reduction using mathematical

programming approaches

Optimization of hepatocyte functionality Reduction of complex chemistry Uncertainty analysis & feasibility evaluation Analysis of alternative solutions

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Uncertainty with Unknown Behavior Uncertainty with Unknown Behavior Alternative solutions that spans the range

• f uncertainty

Alternative solutions that spans the range

• f uncertainty

MILP parametric

• ptimization

MILP parametric

• ptimization

model robustness model robustness solution robustness solution robustness

Determine a Set of Alternative Solutions

Known Behavior

• f Uncertainty

Known Behavior

• f Uncertainty

Robust

• ptimization

method Robust

• ptimization

method A set of solutions represent trade-

• ff between

various objectives A set of solutions represent trade-

• ff between

various objectives

Li and Ierapetritou, Comp. Chem. Eng. in press, 2007 (doi:10.1016/j.compchemeng.2007.03.001).

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Multiparametric MILP (mpMILP) Approach

min ( ) . . {0,1}, 1,...,

T l u j

z c D x s t Ax b E x x j k θ θ θ θ θ = + ≥ + ≥ ≤ ≤ ∈ =

BASIC IDEA ∗ One critical region with one starting point ∗ Complete solution is retrieved with different

starting points (parallelization)

mpMILP problem generalized from scheduling under uncertainty ∗ same integer solution ∗ same parametric objective: z*=f(θ) ∗ same parametric solution (continuous variable): x*= f(θ) In any Critical Region of an mpMILP

Li and Ierapetritou, AIChE Jl. 53, 3183, 2007; Ind. Eng. Chem. Res. 46, 5141, 2007.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Illustrating Example

Ierapetritou MG, Floudas CA, 1998

Determine: Task Sequence Exact Amounts of material processed Given: Raw Materials Required Products Production Recipe Unit Capacity Objective: Maximize Profit S1 S2 S3 S4 mixing reaction purification

Demand and Price Only Demand, Price, Processing Time Uncertainty

2 2 1 2 1 1 2 2

Profit 88.55 49.07 0.25 20 1.2 0.01 θ θ θ θ θ θ = + − + − +

2 1 2 1

Profit 88.55 44 25.16 20 θ θ θ = + − −

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

minimize H or maximize ∑price(s)d(s,n) subject to ∑wv(i,j,n) ≤ 1 st(s,n) = st(s,n-1) – d(s,n) + ∑ρP∑b(i,j,n-1) + ∑ρc∑b(i,j,n) st(s,n) ≤ stmax(s) Vmin(i,j)wv(i,j,n) ≤ b(i,j,n) ≤ Vmax(i,j)wv(i,j,n) ∑d(s,n) ≥ r(s) Tf(i,j,n) = Ts(i,j,n) + α(i,j)wv(i,j,n) + β(i,j)b(i,j,n) Ts(i,j,n+1) ≥ Tf(i,j,n) – U(1-wv(i,j,n)) Ts(i,j,n+1) ≥ Tf(i’,j,n) – U(1-wv(i’,j,n)) Ts(i,j,n+1) ≥ Tf(i’,j’,n) – U(1-wv(i’,j’,n)) Ts(i,j,n) ≤ H, Tf(i,j,n) ≤ H

Duration Constraints Demand Constraints Allocation Constraints Capacity Constraints Material Balances Objective Function

Ierapetritou MG, and Floudas CA, Ind. & Eng. Chem. Res., 37, 11, 4341, 1998

Uncertainty with Known Behavior: Robust Optimization

Scenario-based Robust Stochastic Programming

Requires some statistic knowledge of the input data Optimization of expectations is a practice of questionable validity Problem size will increase exponentially with the number of uncertain parameters

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina ⎣ ⎦

⎣ ⎦

l S m t lt l l m lm S M t S M S t S m m lm

p x a x a x a

l l l l l l l l l l l l

≤ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ Γ − Γ + +

∑ ∑

∈ ∈ Γ = ⊆ ∪

| | ˆ ) ( | | ˆ max

} \ , | | , | } { {

Robust Counterpart Optimization

Soyster’s, Soyster (1973)

Soyster’s Ben-Tal and Nemirovski’s Bertsimas and Sim’s

• Linear
• No flexibility
• Most pessimistic
• Nonlinear
• Flexibility
• Relative smaller number of

variables and constraints

• Linear
• Higher flexibility
• Relative larger number of

variables and constraints

Ben-Tal and Nemirovski’s, Ben-Tal and Nemirovski (2000); Lin, Janak et al. (2004) Bertsimas and Sim’s, Bertsimas and Sim, 2003

2 /

2

, |] | , 1 max[ Pr

Ω −

= ≤ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + ≥ +∑

∑

e p p y b x a

l l k k lk m m lm

κ κ δ

l l M m m lm lm M m m lm

p p u a a x a

l l

ˆ ) ˆ ( − ≤ + + ∑

∑

∈ ∉

Efficient alternative to scenario based robust stochastic programming

] ˆ , ˆ [ ~

lm lm lm lm lm

a a a a a + − ∈

Find solution which copes best with the various realizations of uncertain data

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

700 625 625 625 Continuous variable 1239 1167 1167 1167 constraints 216 216 216 216 Binary variable 254 43.4 4.8 infeasible 150 4.2 CPU time p≤0.5 p≤0.625 p≤0.75 k=75%

• Probability of constraint

violation 939.12 1005.5 1052.50

• 939.12

1052.50

• bjective

Г=1 Г=0.5 Г=0 Bertsimas and Sim Ben-Tal Soyster Deterministic

Comparison for the robust courterpart formulations for processing time uncertainty

• 15% variability for all the processing time
• 72 hours horizon, 24 event points

S1 S2 S3 S4 mixing reaction purification

Illustration

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina ∗Study the effect of the different uncertainties

∗ Provide an efficient way to look up the reactive schedule with the realization of uncertainty (e.g., rush order, machine breakdown)

Parametric and Robust Solution

Parametric Solution Robust Counterpart Solution

∗ Provide an effective way to generate robust preventive schedule with boundary information on uncertainty (e.g., processing time variability)

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Uncertainty in Hepatocyte Functionality

How can we use these techniques to deal with

experimental variability?

In many cases experimental error is more than 100%

How can we analyze the results?

Is the results an artifact of uncertainty?

How can we move beyond experimental error?

Can we determine which parameters are more important and what experiment to do next?

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Aim: Identify the flux distributions for optimal urea production that can fulfill metabolites balances and flux constraints

Unit: µmol/million cells/day

1 min max

: :

urea N ij j j j j j

Max Z v Subject to S v i M v v v j K

=

= = ∀ ∈ ≤ ≤ ∀ ∈

∑

Single-level Optimization: Maximize Urea Secretion

> 2 fold 2.35±0.52 LPAA > 15 fold 0.17±0.24 LIP > 3 fold 1.32±0.69 HPAA > 10 fold 6.81 0.23±0.43 HIP Increase Optimal Value Experimental Data*

*Chan & Yarmush et al (2003) Biotechnol Prog

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Example of Multiple Solution in 2-D

Feasible Region x ≥ 0 y ≥ 0

• 2 x + 2 y ≤ 4

x ≤ 3 Subject to: Minimize x - y

Multiple Optimal Solutions!

4 1

x

3 1 2

y

2 3

1/3 x + y ≤ 4

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Finding all Solutions

Nathan D, Price et al., 2004, Nature Review: Microbiology

A recursive MILP problem that has a set

• f constrains for changing the basis and

identifying a new extreme point A recursive MILP problem that has a set

• f constrains for changing the basis and

identifying a new extreme point

Lee et al., 2000, Computer and Chemical Engineering

T

Z z s t B z q z = = ≥ m in . . α

α

−

∈ ∈ −

= = ≥ ≤ − = − ≤ ≤ ∈ + ≤ ∈ ≥

∑ ∑

K 1 k

K T i i NZ k i i NZ i i K 1 i i

Z z s t Bz q y 1 w NZ 1 k 1 2 K 1 z Uw i I y w 1 i NZ z m in . . , , ,..., , ,

(MILP)

Question: How can you determine all solutions? Question: How can you determine all solutions?

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Flux distributions including glucose production (left) & without glucose production (right)

MILP Model: Application to Hepatocytes

0.327 0.327 0.509 0.509 R(1,72) 0.806 0.806 0.178 0.178 0.806 0.806 0.178 0.178 R(7,6) 3.76 3.76 3.76 3.76 3.76 3.76 3.76 3.76 R6-R7 D8 D7 D6 D5 D4 D3 D2 D1 1 1 7 2

(1, 7 2 ) v R v v = +

7 6 7

) 6 , 7 ( v v v R − =

Enumerate Eight different flux distributions flux distributions that satisfy mass balance and all constraints with the same value of maximal urea production.

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Develop Patient Specific Treatment

ROBUST SOLUTION CONSIDERING VARIABILITY

* All fluxes are in µmol/million cells/day

LINEAR VARIABILITY IN EXTRACELLULAR FLUXES

• All 19 amino acids are indispensable for

maximum function

• Valine and Isoleucine are required at higher

concentrations

1 2 3 4 5 6 7

Thr Lys Phe Gln Met Val Isoleu

Amino Acids Supplementation Concentration (mM)

Point D of Pareto Set Two Stage Approach

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina

Acknowledgements

Kai He Patricia Portillo Eddie Davis Hong Yang Mehmet Orman

Financial Support: NSF, ONR, PRF, ExxonMobil

Zukui Li

Our Web Page: http://sol.rutgers.edu/staff/marianth

Beverly Smith George Saharidis Zhenya Jia Vidya Iyer Yijie Gao

Marianthi Ierapetritou PASI: August 12-21, 2008, Mar del Plata, Argentina