


Full Text  
PAGE 1 MATHEMATICAL MODELING AND OPTIMAL EXPERIMENTAL DESIGN IN SYSTEMS BIOLOGY BY JONATHAN KEEGAN STATZ A Thesis Submitted to the Division of Natural Sciences New College of Florida in partial fulfillment of the requirements for the degree Bachel or of Arts in Applied Mathematics/Biology Under the Sponsorship of Necmettin Yildirim Sarasota, Florida May 2011 PAGE 2 ii Acknowledgements I would like to thank my family members, especially my sister Kaitlin, for their help during the writing process and for reading over many drafts. Thank you to David Mullins and Chris Hart for agreeing to be members of my committee and evaluating my work. Thank you to John Correa, Catherine Martini, Alexander Salisbury, Casey Henderson, and Nicole Robinson for their help i n grammar corrections and editing. Finally, I would like to thank Necmettin Yildirim for the many hours of assistance and support in finishing this thesis. PAGE 3 iii Table of Contents I. Front matter . . . . . . ii 1. Table of Contents . . . . . .iii 2. List of Figures . . . . . .v 3. List of Tables . . . . . vii 4. Abstract . . . . . . viii II. Ma in Text . . . . . . .1 1. Introduction . . . . . . .1 2. Background . . . . . . .3 2.1. Chemical Kinetics . . . . .3 2.2. Reaction Rates and Rate Laws . . . 3 2.3. Enzyme Kinetics . . . . .8 2.4. Enzyme Biology . . . . 8 2.5. The Michae lis Menten Equation . . . .11 2.6. Enzyme Inhibition . . . . .17 2.7. Hill Function and Sigmoidal Kinetics . . .26 2.8. Simulation of Reaction Kinetics . . . .27 2.9. More Complex Reaction Networks . . 30 3. Methods . . . . . . .39 3.1. Model Discrimination . . . . .39 3.2. Symmetric Matrices and Quadratic Forms . . .41 3.3. Observability and the Observability Gramian . 44 3.4. Multiple Models . . . . 49 PAGE 4 iv 3.5. Derivation of P from L 2 Norm . . . 50 3.6. Optimum Initial Conditions . . . 51 3.7. Application of Method to Enzyme Inhibition Models . 53 4. Results . . . . . . .60 4.1. Differentiation of Two Feedback Models from Experimental Data . . . . . . 60 4.2. Generation of Artificial Experimental Da ta . . 63 4.3. Optimal Experimental Design . . . .67 4.4. Complete Experimental Protocol . . . 72 5. Discussion and Conclusion . . . . .74 III. Back Matter . . . . . . .78 1. Bibliography . . . . . .78 PAGE 5 v List of Figures 1. Product Concentration vs. Rate Plot for Zero Order Reaction . .5 2. Reactant Concentration vs. Rate Plot for First Order Reaction . 6 3. Reactant Concentration vs. Rate Plot for Second Order Reaction . .7 4. Energy Levels of Catalyzed and Uncatalyzed Reactions . . .10 5. Concentrations Over Time for Michaelis Menten Kinetics . . 16 6. Competitive Inhibition Lineweaver Burk Plot . . . .21 7. Uncompetitive Inhibition Lineweaver Burk Plot . . . .23 8. Non competiti ve Inhibition Lineweaver Burk Plot . . 24 9. Mixed Inhibition Lineweaver Burk Plot . . . .25 10. Michaelis Menten vs. Hill Function Comparison . . . .27 11. Time vs. Product Concentration for Zero Order Reaction . . .28 12. Enzyme System Concentrations Over Time . . . .29 13. Mutual Activati on and Mutual Inhibition Bistability Curves . .32 14. Reaction Network for Mutual Inhibition . . . 33 15. Mutual Inhibition Over Time with Bistability . . . 34 16. Reaction Network for Damped and Sustained Oscillations . . 35 17. Example of Damped Oscillations . . . . .36 18. Example o f Sustained Oscillations . . . . .37 19. Difference Between Outputs of Competitive and Uncompetitive Inhibition Models Over Time . . . . . . 58 20. Mechanisms for Models 1 and 2 . . . . .60 21. Time Series for Models 1 and 2 for Initial Values of Zero with Data. . .66 PAGE 6 vi 22. Time Series f or Models 1 and 2 for Optimal Initial Values with Data. . 70 23. Difference Between the Outputs of Models 1 and 2 . . .71 PAGE 7 vii List of Tables 1. Rate Enhancement by Selected Enzymes . . . .9 2. Parameter Values for Optimal Initial Condition Design for Inhibition Models . . . . . . . . 57 3. Parameter Values for Optimal Initial Condition Design for Models 1 and 2 .63 4. Artificially Created Data Points From Initial Values of Zero . .64 PAGE 8 viii MATHEMATICAL MODELING AND OPTIMAL EXPERIMENTAL DESIGN IN SYSTEMS B IOLOGY Jonathan Keegan Statz New College of Florida, 2011 ABSTRACT Mathematical modeling can sometimes be used to better understand how complex reaction networks behave in biological systems. However, a problem often encountered is the existence o f multiple models developed from different interaction mechanisms that fit experimental data equally well. An optimal experiment such that these mechanisms can be discriminated is critical in order to determine the correct mechanism. Mathematical models can be used in designing these experiments that will invalidate the incorrect mechanisms and reduce the time and cost of using multiple, non optimal experiments to achieve the same goal. In this study a discrimination method of initial condition optimizat ion to maximize the differences between model outputs is described and detailed. Then it is applied to two inhibition models in enzyme kinetics to show how the method works. It is then used on a hypothetical system with two different types of positive fe edback loops to develop the complete experimental protocol for designing an optimized experiment directly from experimental data. This method shows that only a total of two time series data sets are needed to differentiate two competing models. Necmettin Yildirim Division of Natural Sciences PAGE 9 1 1. Introduction Studies in the fields of biology look at reactions at every level of organization, from molecular to cellular to organismic. The traditional methods of studying these reactions at the cellular level t ypically involve a single perturbation and a single measurement of how this perturbation affects an isolated system. These traditional methods include reductive methods, those that simplify reactions and systems down to their fundamental parts. The appro ach in systems biology is to study perturbations at the system wide level using integrative or holistic methods. These methods study how perturbations affect other systems, and every component is considered as an interconnected, interacting part of one wh ole system. In both methods mathematical models can be used to describe the mechanisms that express the interactions of these systems. Mathematical modeling is widely used to better understand the complex interaction mechanisms in systems biology. With models the behavior and dynamics of complex systems can be better predicted and investigated. Mathematical models can help researchers design better experiments for the analysis of a reaction mechanism, and analysis of the behavior of simple networks can allow for predictions about the behavior of similar, more complex networks to be made. When using mathematical models to investigate the complex interaction mechanisms of a system there are many instances where there can be two or more competing models de scribing structurally different mechanisms, all of which are of similar complexity and fit a given set of experimental data equally well. Without further PAGE 10 2 methods of optimization, it would take multiple runs of many different experiments to properly identi fy the correct model, a process both costly and time consuming. Therefore, experiments in these situations should be designed such that they provide the most efficient means of discriminating between the various mechanisms involved. An optimized experime nt also helps in reducing the overall cost and time put into a study, as this new experiment should provide the information that standard methods would require multiple experiments to produce. The purpose of this study is to focus on the method of model d iscrimination via optimal initial condition design. The discrimination method will be applied to two different types of enzyme inhibition models (competitive and uncompetitive inhibition) in enzyme kinetics as a proof of concept. It will then be applied t o two competing positive feedback mechanisms to develop an experimental protocol starting from experimental data. The outline of this thesis is as follows: In Chapter 2 the details required to understand the workings of chemical kinetics, rate laws, and enzymes, inhibitors, and their kinetics will be given. In Chapter 3 the method of optimal experimental design, as well as required background information and theory, will be explained in full. In Chapter 4 the results of this method will be demonstrated for two feedback systems and the complete experimental protocol will be described. PAGE 11 3 2. Background This section will include the necessary background, mathematical and biological, that will be required to understand the methods that will be described and used. Included will be a background on reaction kinetics, enzyme kinetics, and the Mic haelis Menten equation. There will also be detailed sections on more complex reaction networks and their simulations. 2.1. Chemical Kinetics Chemical kinetics is a branch of physical chemistry that studies the dependence of chemical reaction rates on reac tant concentration. Reactions can be classified by one of two characteristics: molecularity or order. Molecularity refers to how many molecules are altered in the reaction, or how many molecules are in the activated complex. Molecularity is always an i nteger and seldom greater than three, meaning almost all reactions are unimolecular, bimolecular, or termolecular. The order of a reaction mathematically describes a reaction's kinetics and defines the number of concentration terms multiplied to get a rat e law for the reaction. These rate laws can be used to determine how long a reaction will take and what the final yield of the reaction will be (Leskovac, 2003). 2.2. Reaction Rates and Rate Laws The individual steps of larger reactions tend to have sma ller orders. As per the law of mass action, the rates of these smaller reactions, and all following rate laws, are directly proportional to the product of the concentrations of each reactant. Two general mechanisms and rate equations used to illustrate t his law are: PAGE 12 4 For single molecules of each reactant: A + B k P d P [ ] dt = k A [ ] B [ ] For reactions where there are m molecules of A and n molecules of B : m A + nB k P d P [ ] dt = k A [ ] m B [ ] n where brackets indicat e the concentration of a specific component, k is a rate constant, and t is time. In zero order reactions the rate of the reaction is constant and unaffected by reaction concentration. These reactions can be demonstrated with k # P with the rate law v = d P [ ] dt = k ( 2 1 ) where v is the production rate of P [ ] P [ ] is the product concentration at time t and k is the rate constant. As the mechanism shows, a zero order reaction indicates the presence of a catalyst and very large amounts of reactant, such that the for mation of the product has little to no effect on the amount of reactant and the reaction continues at a constant rate unimpeded until the experiment is stopped. PAGE 13 5 Figure 1. This figure shows the concentration of product vs. rate reaction for a zero order reaction described by Eq. ( 2 1 ) For a first order irreversible reaction such as A k P the rate law is known to be v = d P [ ] dt = d A [ ] dt = k A [ ] = k A 0 [ ] P [ ] ( ) ( 2 2 ) where A [ ] is the reactant concentration at time t A 0 [ ] is the total concentration of A in all its forms (here that is A and P ), and k is the first order rate constant with units 1 /time PAGE 14 6 Figure 2. A graphical representation of the concentration of reactant vs. rate reac tion for a first order reaction as described by Eq. ( 2 2 ) The rate law for a second order irreversible reaction, such as A + B k P is v = d P [ ] dt = d A [ ] dt = k A [ ] B [ ] ( 2 3 ) where k is the second order rate constant with unit concentration/time Th rough the use of the mass conservation equations A 0 [ ] = A [ ] + P [ ] and B 0 [ ] = B [ ] + P [ ] Eq. ( 2 3 ) can be written as v = k A 0 [ ] P [ ] ( ) B 0 [ ] P [ ] ( ) where A 0 [ ] and B 0 [ ] are the total concentrations of A and B in all their forms, respectively. PAGE 15 7 Figure 3. This figure displays the concentration of reactant vs. rate reaction for a second order reaction as described by Eq. ( 2 3 ) where B [ ] is a second molecule of A [ ] as per the methods of the law of mass action. When the reactant B is a second molecule of A the second order reaction can be written as 2 A k P The rate law for this reaction can be written as v = k A [ ] 2 = k A 0 [ ] P [ ] ( ) 2 However, most reactions are reversible, and the rate equation must be altered to in clude the reverse reaction. A reaction involving a single reactant can be represented with A k 1 k 2 # P The rate law of the production of P [ ] is the difference of the forward and backward rates: PAGE 16 8 v = d P [ ] dt = d A [ ] dt = k 1 A [ ] k 2 P [ ] = k 1 A 0 [ ] k 1 + k 2 ( ) P [ ] ( 2 4 ) where k 1 and k 2 are the forward and reverse rate constants, respectively (Leskovac, 2003). This can also be written in terms of the rate of change of A [ ] as d A [ ] dt = k 1 + k 2 ( ) A [ ] k 2 A 0 [ ] 2.3. Enzyme Kinetics Enzyme kinetics is the type of chemical kinetics that focuses on enzyme catalyzed biochemical reactions. Th ese reactions tend to not have simple orders, but the individual steps in a reaction do, allowing for methods of analysis similar to those used for chemical kinetics. Zero order reactions can occur in enzyme kinetics when there is a very large amount of s ubstrate present (Bisswanger, 2008), but most of the reactions in this field tend to have very complicated rate equations. Biochemical reactions are often modeled as a set of ordinary differential equations (ODEs) with or without nonlinearities that descr ibe the reactions involving a specific component in a stoichiometric setting, usually in physiological conditions where the reactions are in some sort of equilibrium. Using quasi steady state approximations the number of dynamic variables and parameters c an be reduced to achieve a simpler model, as will be explained later. A system of ODEs can be used to describe reactant concentrations over time, and from simpler models complex models involving multiple reactants and products in a single network can be d escribed. 2.4. Enzyme Biology Enzymes are catalytic proteins that increase the rate at which a reaction occurs by lowering the activation energy required. This is achieved when the enzyme binds with PAGE 17 9 the substrate(s) to form an enzyme substrate complex. Table 1 shows this increase in reaction rate for a selection of enzymes. They are highly specific in both the reactions they catalyze and the substrates to which they bind, usually only catalyzing a single reaction or a set of closely related reactions. This portion, where the catalytic residues reside, the substrate binds, and the reaction occurs, is known as the active site. They may also have smaller binding sites for smaller substrates, products of the catalyzed reaction, or other effector molecules whose binding can increase or decrease the enzymatic activity. Table 1. Rate enhancement by some enzymes. Enzyme Uncatalyzed Rate ( k un s 1 ) Catalyzed Rate ( k cat s 1 ) Rate Enhancement ( k cat s 1 / k un s 1 ) OMP decarboxylase 2.8 10 16 39 1.4 10 17 Staphylococcal nuclease 1.7 10 13 95 5.6 10 14 AMP nucleosidase 1.0 10 11 60 6.0 10 12 Carboxypeptidase A 3.0 10 9 578 1.9 10 11 Ketosteroid isomerase 1.7 10 7 66000 3.9 10 11 Triose phosphate isomerase 4.3 10 6 4300 1.0 10 9 Chorismate mutase 2.6 10 5 50 1.9 10 6 Carbonic anhydrase 1.3 10 1 1 10 6 7.7 10 6 Abbreviations: OMP, orotidine monophosphate; AMP, adenosine monophosphate. ( Berg, 2006). Enzymes work by lowering the activation energy of a reaction, G A P the amount of energy required for the conversion of reactants to products to begin. They do not alter the equilibrium concentrations of a chemical reaction significantly, as the equilibrium is based upon the change i n Gibbs standard free energy, G O PAGE 18 10 Figure 4. A plot comparing free energy during the reaction for both the catalyzed and uncatalyzed reaction. The difference in the required activation energy can be clearly seen. At equilibri um the equilibrium constant is known to be G O = # R T l n K M where R is the gas constant 8.315 J / m ol K T is the temperature at which the reaction occurs in Kelvin, and K M is the equilibrium constant for the equation, which will be derived further below (Berg, 2006). The basic mechanism for an enzyme to lower the activation energy is to bind to the substrate to form an enzyme substrate complex. This complex makes the conditions in which the reaction occurs m ore favorable for the formation of the product. This improvement in conditions occurs in one of several ways; examples include: the enzyme can cause a conformational change in the substrate, making it closer to the transition state form and reducing the energy needed to reach the transition state; it can reduce the entropy of the reaction by bringing the substrates together into the correct conformation, reducing the energy needed to orient them correctly; or it can increase the temperature at which the r eaction occurs, speeding up the reaction (although temperature must be regulated to prevent denaturing the enzyme). Once the product is formed it dissociates PAGE 19 11 from the complex and the enzyme moves on to catalyze another reaction. The mechanism for this is shown here: E + S E S E + P where E is the enzyme, S is the substrate, ES is the enzyme substrate complex, and P is the product formed in the reaction. There are several factors that can affect the rate of enzymatic reactions; for example, the y are heavily dependent on environmental factors such as temperature and pH. At high enough temperature, enzymes, like all proteins, can become denatured, losing their shape and activity. A similar result can occur from exposure to extreme pH levels in e ither direction, or to very high salt concentrations. Depending on the enzyme, denaturation may be irreversible. Other molecules can also affect the activity of an enzymatic reaction; these molecules can include inhibitors, the reactions own products and substrates, and allosteric effectors (Berg, 2006). 2.5. The Michaelis Menten Equation The most commonly used mathematical model for biochemical reaction kinetics is the Michaelis Menten equation, which relates the concentrations of a substrate to the reac tion rates in order to describe the rates of irreversible enzymatic reactions. A common Michaelis Menten equation is derived from the following reaction mechanism: E + S k 1 k 2 # E S k 3 P + E The rate determining step of this reaction depends on the speci fic enzymes and substrates involved; depending on their identities the slowest step could involve: a conformational change in the enzyme or substrate during binding; the chemical reaction needed to create the product; or possibly a conformational change i n the enzyme in order to release the PAGE 20 12 newly formed product. In general it is the step that requires the greatest activation energy to begin or that with the transition state with the highest free energy. Whichever step is the slowest, it will determine th e overall kinetics of the catalyzed reaction (Bisswanger, 2008). According to mass action kinetics, the system of differential equations describing the changes in concentrations of each component over time is d S # $ dt = k 2 E S # $ % k 1 S # $ E # $ d E S # $ dt = k 1 S # $ E # $ % k 2 + k 3 ( ) E S # $ d E # $ dt = k 2 + k 3 ( ) E S # $ % k 1 S # $ E # $ d P # $ dt = k 3 E S # $ The sum of the rates of change for the enzyme in all its forms is equal to zero, that is, d E # $ dt + d E S # $ dt = 0 Integration of this equation provides the equation for the conservation o f the total enzyme concentration E 0 # $ as E 0 # $ = E # $ + E S # $ ( 2 5 ) Similar to the enzyme concentration, the sum of the rates of change for the substrate forms is d S # $ dt + d E S # $ dt + d P # $ dt = 0 PAGE 21 13 which provides the mass conservation equation for the total substrate S 0 # $ as: S 0 # $ = S # $ + E S # $ + P # $ ( 2 6 ) Through the use of Eqs. ( 2 5 ) and ( 2 6 ) the system of differential equations for the rates can be reduced from four to two: d E S # $ dt = k 1 S # $ E # $ % k 2 + k 3 ( ) E S # $ ( 2 7 ) d P # $ dt = k 3 E S # $ ( 2 8 ) One form of the Michaelis Menten equation is based on the two assumptions: (1) k 3 k 2 which simplifies Eq. ( 2 7 ) to d E S # $ dt = k 1 S # $ E # $ % k 2 E S # $ ( 2 9 ) and (2) the enzyme substrate complex, free substrate, and free enzyme are in thermodynamic equilibrium, which simplifies Eq. ( 2 9 ) to d E S # $ dt = 0 = k 1 S # $ E # $ % k 2 E S # $ which results in K M = k 2 k 1 = E # $ S # $ E S # $ ( 2 10 ) Solving Eq. ( 2 5 ) for E [ ] and substituting the result into Eq. ( 2 10 ) allows for the concentration of the enzyme substrate complex to be written as E S # $ = S # $ E 0 # $ K M + S # $ ( 2 11 ) PAGE 22 14 Substituting Eq. ( 2 11 ) into Eq. ( 2 8 ) gives d P # $ dt = v = V m a x S # $ K M + S # $ ( 2 12 ) where V m a x is the maximum rate of the reaction: V m a x = k 3 E 0 # $ ( 2 13 ) For enzymatic reactions, it is assumed E 0 # $ S 0 # $ ; this a ssumption implies that E S # $ S # $ + P # $ which reduces Eq. ( 2 6 ) to S # $ = S 0 # $ % P # $ ( 2 14 ) The differential equation in Eq. ( 2 12 ) with Eq ( 2 14 ) substituted for S # $ allows for an equation describing the concentration of the product over time in terms of only the product However, there are some problems with this method when doing experiments in living systems: within real life examples of biochemical and cellular systems the conditions may differ significantly from the ideal conditions upon which equilibrium kinetics r elies. An assumption made is that all the rate constants are fixed and that the k 3 reaction must be much lower than the reaction containing both k 1 and k 2 However, in most biological systems the rate constants actually vary within a large value range du ring reactions (Berg, 2006). In order to derive a better equation describing the rate of the reaction another assumption must be made: the quasi steady state assumption This will allow for the derivation of the true Michaelis Menten equation. The quas i steady state assumption assumes that the concentration of the enzyme complex changes more slowly than those of the substrate and the product and is nearly PAGE 23 15 constant throughout the reaction; in situations where this assumption is true many, but not all, of the parameters are held constant. This allows the rate of change for the concentration of the enzyme substrate concentration, Eq. ( 2 7 ) to be set to zero; that is, d E S # $ dt = 0 = k 1 S # $ E # $ % k 2 + k 3 ( ) E S # $ ( 2 15 ) The Michaelis Menten constant taken from Eq. ( 2 15 ) is now K M = k 2 + k 3 k 1 = E # $ S # $ E S # $ ( 2 16 ) As seen in Eq. ( 2 16 ) the Michaelis constant is increased by k 3 / k 1 compared to the constant in Eq. ( 2 10 ) Solving Eq. ( 2 5 ) for E [ ] and substituting the result into Eq. ( 2 15 ) allows for E S [ ] in terms of substrate concentration and initial enzyme concentrations to be solved for: E S # $ = k 1 S # $ E 0 # $ k 1 S # $ + k 2 + k 3 When this is substituted into Eq. ( 2 8 ) the following equation is obtained d P # $ dt = v = k 3 E S # $ = k 3 S # $ E 0 # $ S # $ + k 2 + k 3 k 1 This can be simplified to the Michaelis Menten equation: v = V m a x S # $ K M + S # $ ( 2 17 ) PAGE 24 16 Here V max is the same as Eq. ( 2 13 ) (Bisswanger, 2008). For both cases the plot of the rate equations as a function of S # $ gives an increasing hyperbolic curve that approaches but never reaches V max Figure 5. This figure s hows the concentratio ns of the substrate, product, enzyme, and enzyme substrate complex during the reaction time under the quasi steady state assumption made by the Michaelis Menten method. It is also assumed that E 0 [ ] S 0 [ ] When a reaction has just begun there is no product in the system and the concentration of the substrate is equal to the total concentration of the substrate from the law of mass action; that is, S [ ] S 0 [ ] This means that for a small time frame near the beginning of the reaction the Michaelis Menten equation can be written as v = V m a x S 0 [ ] K M + S 0 [ ] ( 2 18 ) This allows for the study of reactions that occur very rapidly. Plotting the inverse of Eq. ( 2 18 ) : 1 v = K M V m a x S 0 [ ] + 1 V m a x ( 2 19 ) S ES E P Reaction Time Concentration PAGE 25 17 will give the Lineweaver Burk plot, which allows for the values of the combined rate parameter and maximum rate to be easily t aken from a plot of 1 / S 0 [ ] versus 1/ v 2.6. Enzyme Inhibition This section will discuss the types of inhibitors and their interactions with the enzyme; specific inhibitor kinetics will be discussed later. There are four types of rever sible inhibitors, which can be differentiated from each other through analysis of their kinetics as more inhibitor is introduced into the system. In competitive inhibition the inhibitor is similar in shape to the substrate and binds to the empty active sit e of the enzyme, competing with the substrate to form a complex first. The effects of this inhibitor can be overcome with very high concentrations of the substrate. In uncompetitive inhibition the inhibitor binds to the enzyme after it has already been b ound to the substrate, delaying the release of the substrate or product. The ability of this type of inhibitor to bind to the complex may be dependent on a conformational change that only occurs when the substrate and enzyme bind, forming a new binding si te for the inhibitor. Another inhibition type is mixed inhibition. Here the inhibitor does not interfere with substrate binding, but the binding of the inhibitor or substrate affects the ability of the other to bind to an enzyme, via, for example, throug h the inducing of a conformational change in the enzyme's active site. Although it is possible for mixed inhibition to occur at the active site, this type of inhibition is more likely to be due to allosteric regulation; that is, binding occurs elsewhere o n the enzyme. The last type of reversible inhibition is non competitive inhibition. Non competitive inhibition is a special case of mixed inhibition, where the rate at which the inhibitor binds to the enzyme is equal to the rate at which it binds to the PAGE 26 18 enzyme substrate complex. Here an inhibitor can bind to the enzyme, though not to the active site, but does not interfere with the binding of the substrate. The amount of inhibition that occurs depends entirely on the inhibiter concentration. If, after the inhibitor has bound to the enzyme, there is still some residual enzymatic activity in the enzyme substrate inhibitor complex, resulting in a tiny amount of product formation, the inhibition effect is known as partial inhibition (Berg, 2006). When data from an experiment involving inhibition is collected, normal methods to determine the type of inhibition include multiple runs with different concentrations of the inhibitor, such that the behavior of the system can be clearly matched to a known inhibition behavior. This would involve a large amount of work, time, and resources as each inhibition concentration, of which there must be at least two, also requires multiple runs with many initial substrate concentrations to correctly identify the type of inhib ition occurring. This method can be demonstrated using the enzyme inhibitor types described earlier and their corresponding Michaelis Menten equations and Lineweaver Burk plots. As will be seen, each inhibition type has its own unique behavior of the Line weaver Burk plot. From earlier it is known that the mechanism and the rate equation for product formation in an uninhibited enzyme reaction are E + S E S E + P and v = V m a x S [ ] K M + S [ ] For competitive inhibition the mechanism is PAGE 27 19 k 1 k 3 E + S + I ES E + P EI k 2 k i k i w here I and EI are the inhibitor and inhibitor enzyme complex, respectively. The total concentration of the enzyme in all its forms is assumed to be constant, represented by: E 0 [ ] = E [ ] + E S [ ] + E I [ ] ( 2 20 ) The rate equations for each of these forms of the enzyme can be written as d E [ ] dt = k 1 E [ ] S [ ] + k 2 E S [ ] + k 3 E S [ ] k i E [ ] I [ ] + k i E I [ ] ( 2 21 ) d E S [ ] dt = k 1 E [ ] S [ ] k 2 + k 3 ( ) E S [ ] ( 2 22 ) d E I [ ] dt = k i E [ ] I [ ] k i E I [ ] ( 2 23 ) where I [ ] is the inhibitor concentration and E I [ ] is the concentration of the enzyme inhibitor complex. The rate equation for the product can be written as d P [ ] dt = k 3 E S [ ] ( 2 24 ) Assuming d E I [ ] dt = 0 that is, it is at equilibrium, allow s for the inhibitor dissociation constant K i for binding to the enzyme to be solved for as K i = k i k i = E [ ] I [ ] E I [ ] PAGE 28 20 Setting Eq. ( 2 22 ) equal to zero will provide Eq. ( 2 16 ) and solving for E [ ] gives E [ ] = k 2 + k 3 ( ) E S [ ] k 1 S [ ] = K M E S [ ] S [ ] ( 2 25 ) Substituting Eq. ( 2 25 ) into Eq. ( 2 23 ) assuming Eq. ( 2 23 ) is in equilibrium, and solving for E I [ ] provides E I [ ] = K M I [ ] E S [ ] k i k i S [ ] = K M I [ ] E S [ ] K i S [ ] ( 2 26 ) Substituting Eqs. ( 2 25 ) and ( 2 26 ) into Eq. ( 2 20 ) and solving for E S [ ] gives E S [ ] = E 0 [ ] S [ ] S [ ] + K M 1 + I [ ] K i # $ % & ( 2 27 ) Inserting Eq. ( 2 27 ) into Eq. ( 2 24 ) allows for the rate of production of P to be found: v = d P [ ] dt = V m a x S [ ] S [ ] + K M where = 1 + I [ ] K i # $ % & and V m a x = k 3 E 0 [ ] As can be seen, this type of inhibition increases the apparent K M value seen in the uncatalyzed reaction (Eq. ( 2 18 ) ) without affectin g the maximum velocity of the reaction. A comparison of the Lineweaver Burk plots of uninhibited and competitively inhibited systems would show that the slope of the inhibited mechanism becomes steeper as the inhibitor concentration is increased, pivoting about the y axis intercept that defines the PAGE 29 21 E + S ES + I E + P k 1 k 2 k 3 ESI K i V max value. Fig. 6 shows a comparison of the Lineweaver Burk plot for the uninhibited enzyme (Eq. ( 2 19 ) ) with the Lineweaver Burk plot of competitive inhibition: 1 v = K M V m a x S 0 [ ] + 1 V m a x Figure 6. This figure shows the effect of a competitive inhibitor on the Lineweaver Burk plot and the values of the constants compared to an uninhibited enzyme. The dashed line shows the plot for the uninhibited enzyme. For the second reversible inhibition type, uncompetitive inhibition, the mechanism is where ESI is the enzyme substrate inhibitor complex. The new rate equation for product formation is PAGE 30 22 v = V m a x S [ ] K M + S [ ] where = 1 + I [ ] K i where K i is the inhibitor dissociation constant for binding to the enzyme substrate complex: K i = E S [ ] I [ ] E SI [ ] As can be seen, this type of inhibition decreases both the apparent K M value and the maximum velocity of the reaction. A comparison of the Lineweaver Burk plots of uninhibited and uncompetitively inhibited systems would show that the slope of the inhibited mechanism is unchanged as the inhibitor concentration increased. Instead, a parallel plot with axis intercepts higher in value as inhibitor concentration increases would be present. It is important to note that as I [ ] increases, K M and V max are decreased by an equal amount compared to each other. Fig. 7 shows a comparison of the Lineweaver Burk plot for the uninhibited enzyme (Eq. ( 2 19 ) ) with the Lineweaver Burk plot of uncompetitive inhibition: 1 v = K M V m a x S 0 [ ] + V m a x PAGE 31 23 E+S + I ES + I E+P EI+S k 1 k 2 k 3 K i ESI K SS K i Figure 7. This figure displays the effect of an uncompetitive inhibitor on the Lineweaver Burk plot and the values of the constants compared to an uninhibited enzyme. The dashed line shows the plot for the uninhibited enzyme. Non competitive in hibition can be represented with Here the new rate equation for product formation is v = V m a x S [ ] K M + S [ ] As will be seen below, non competitive inhibition is a case of mixed inhibition where = This type of inhibitio n decreases the apparent maximum velocity without affecting the K M value of the reaction. A comparison of the Lineweaver Burk plots of uninhibited and non competitively inhibited systems would show that the slope of the inhibited mechanism becomes steeper as the inhibitor concentration increased, pivoting about the x axis intercept that defines the K M value. Fig. 8 shows a comparison of the PAGE 32 24 E+S + I ES + I E+P EI+S k 1 k 2 k 3 K i ESI K SS K i EI+P Lineweaver Burk plot for the uninhibited enzyme (Eq. ( 2 19 ) ) with the Lineweaver Burk plot of noncompetitive inhibition: 1 v = K M V m a x S 0 [ ] + V m a x Figure 8. A graphical representation of the effect of a non competitive inhibitor on the Lineweaver Burk plot and the values of the con stants compared to an uninhibited enzyme. The dashed line shows the plot for the uninhibited enzyme. The last reversible inhibition type is mixed inhibition, which can be shown with: The new rate equation for product formation is v = V m a x S [ ] K M + S [ ] PAGE 33 25 As can be seen this type of inhibition decreases the apparent V max while changing the apparent K M value of the reaction. Which direction the change in the apparent K M occurs depends on whether > or < ; the first would lead to an increasing value, the second a decreasing value. A comparison of the Lineweaver Burk plots of uninhibited and mixed inhibition systems would show that the slope of the inhibited mechanism becomes steeper as the inhibitor concentration increased, pivoting about a point between the x and y intercepts as K M increases (thus moving the x intercept towards zero) and V max decreases (increasing the y intercept) (Berg, 2006). Fig. 9 shows a comparison of the Lineweaver Burk plot for the uninhibited enzyme (Eq. ( 2 19 ) ) with the Lineweaver Burk plot of mixed inhibition: 1 v = K M V m a x S 0 [ ] + V m a x Figure 9. This figure shows the effect of a mix ed inhibitor on the Lineweaver Burk plot and the values of the constants compared to an uninhibited enzyme. The dashed line shows the plot for the uninhibited enzyme. Here > As can be seen for the case of inhibitors, when gi ven experimental data for two or more inhibited reactions, as well as the necessary data to form a model for an uninhibited reaction, it would be possible to not only differentiate between the two inhibited PAGE 34 26 reactions, but to glean some information about th eir equations and mechanisms from given data. However, the amount of data needed would require a large amount of experiments to collect. By studying the kinetics of these enzymatic reactions, large amounts of data can be revealed about an enzyme: the sp ecifics on how the enzyme works and its role in metabolism and homeostasis, how an enzyme's activity is controlled via external signals, and how synthetic signals and molecules can affect its activity. This sort of data is invaluable in medical research, leading to new and more effective drugs and treatments. 2.7. Hill Function and Sigmoidal Kinetics Another important enzyme kinetic model is the Hill Function, which allows for the rate equation of a reaction involving enzymes with multiple substrate bindi ng sites to be defined. This type of reaction can be represented by E + nS E S n E + P Assuming the first step of the reaction is at equilibrium, the equilibrium equation is K d = E [ ] S [ ] n E S n [ ] and the rate equation for production of the product, known as the Hill equation, is d P [ ] dt = v = V m a x S [ ] n K d + S [ ] n ( 2 28 ) In the Hill eq uation the constant K d is the apparent dissociation constant. The parameter n is the Hill coefficient, whose value describes the cooperativity of the substrate binding. A value greater than one indicates a positively cooperative reaction, where the bindi ng of one substrate increases the enzyme's affinity for other substrate molecules, which provides sigmoidal behavior. An n value less than one indicates the opposite: a PAGE 35 27 negatively cooperative reaction where the binding of a substrate lowers the binding a ffinity. When n equals one the reaction is noncooperative: substrate binding has no affect on the enzyme's affinity for further substrate binding and Eq. ( 2 28 ) becomes t he Michaelis Menten equation. A plot of the Hill equation with an n value greater than one, where more than one substrate can bind to an enzyme, is an increasing sigmoidal curve that approaches but does not reach the value of V max A comparison of a Mich aelis Menten curve and a Hill equation curve is shown in the Fig. 10 below. If these reactions are run longer they eventually reach the same value, the shared V max (Leskovac, 2003). Figure 10. This figure depicts a comparison of the plots for general cases of the Michaelis Menten equation and the Hill equation for n >1 from Eqs. ( 2 18 ) and ( 2 28 ) 2.8. Simulation of Reaction Kinetics The ultimate goal of mathematically modeling a system is to study the change in the concentrations of each component in a reaction network over time. This can be achieved by solving the differential equati on(s) once the parameters, initial concentrations, and total concentrations have been fixed. Finding this solution is possible via basic integration methods for simpler models, such as that for a zero order reaction. But for more complex systems, the ana lytical solution is often not possible due to nonlinearities in the system. Therefore, numerical methods must be used to solve these PAGE 36 28 equations. Computer programs can be used for these methods; these programs include Matlab and Maple, which are used for t his study. Via simple integration of the rate equation for the zero order reaction, Eq. ( 2 1 ) with an initial product of P 0 ( ) = 0 the equation describ ing the product concentration at time t is P t ( ) = k t ( 2 29 ) Figure 11. A simulation of the concentration of P over time for k values from 0 to 5 for a zero order reaction from Eq. ( 2 29 ) with P 0 ( ) = 0 For the enzym e kinetics system described previously, E + S k 1 k 2 # E S k 3 E + P ( 2 30 ) the system of equatio ns describing the dynamics of this system is PAGE 37 29 d S [ ] dt = k 2 E S [ ] k 1 S [ ] E [ ] d E [ ] dt = k 2 + k 3 ( ) E S [ ] k 1 S [ ] E [ ] d E S [ ] dt = k 1 S [ ] E [ ] k 2 + k 3 ( ) E S [ ] d P [ ] dt = k 3 E S [ ] The concentrations of each of the four variables over time can be found by solving this system of ODEs numerically after fixing initial conditions for each variable and assigning a va lue to each parameter. Fig. 12 shows the simulation result of this system with initial conditions S 0 [ ] = 2 E 0 [ ] = 2 E S 0 [ ] = 0 and P 0 [ ] = 0 and rate parameters k 1 = 0.5 k 2 = 1 and k 3 = 1 Figure 12. This figure shows the concentrations over time for the four variables of the enzyme kinetic equations. See text for the initial concentrations and rate parameters. I nitially, there has to be non zero initial conditions for both the enzyme and the substrate in order for the reaction ( 2 30 ) to progress, while both the enzyme substrate co mplex and PAGE 38 30 product initial concentrations can begin at zero when t =0. As the reaction progresses the concentration of the substrate S decreases as it is converted to the product P while the product concentration increases in a nearly equal behavior. The concentration of the enzyme E decreases as it binds with the substrate to become the complex EP then returns back to its initial concentration, while the complex concentration increases as it is formed then returns to zero as the substrate is fully conver ted. 2.9. More Complex Reaction Networks Combining, chaining, and linking several of these types of simpler reactions can be used to build more complex systems and lead to far more complex dynamics and long term behaviors. Feedback systems are often used as a means to control the final product concentration at certain substrate levels. In order to prevent problems in the cellular system involving multiple variables the reactions must be carefully regulated to maintain the proper amount of products and rea ctants. This regulation can occur through feedback loops and occasionally feed forward loops. Feedback loops occur when some of the product is fed back into the reaction as a means to control the rate of its production. If the production of the product in some way increases the reaction that created it, by, for example, increasing the substrate binding rate, or by making more substrate available to the enzyme during an equal amount of time, the feedback is positive If the creation of the product in the reaction reduces the rate of product formation the feedback is negative In cellular systems negative feedback is generally used to maintain stability and homeostasis despite outside influences. Some enzyme inhibitors are examples of this type of feedba ck. Positive feedback is used to amplify a small change to cause even larger changes to PAGE 39 31 occur. These sorts of loops occur in morphogenesis, organ development, and blood clotting, as well as many other systems, and are usually controlled by a negative fee dback system, which prevents the changes from getting too large and eventually breaks or overcomes the positive loops (Mitrophanov, 2008). Feed forward loops occur when a system reacts to a change in the external environment without having to be influence d by it first. A system with this kind of regulation is anticipatory, and must react to the change before a feedback system can react to any deviations in the behavior of the system caused by the change in the external environment. Examples of physiologi cal systems that use feed forward loops are gene regulation and transcription, as well as cell growth and division ( Csiksz Nagy, 2009 ). Positive Feedback and Bistability Positive feedback includes two discontinuous types of regulatory systems: (1) mutual activation and (2) mutual inhibition. In mutual activation the synthesis of one reaction component activates a second component, which leads to the increased synthesis of the first component. Mutual inhibition occurs when the first component inhibits th e second component, which in turn promotes the degradation of the first component. In both of these cases a discontinuous switch is created; mutual activation creates a one way switch and mutual inhibition creates a two way switch. For mutual activation a nd the one way switch, when the value of a specific parameter is increased beyond the critical value, the system no longer displays bistability and the final value of the variable being measured will change suddenly and irreversibly to the upper steady sta te; returning the parameter to its previous values will not return the product concentration to its previous value (Fig. 13(a)). Mutual inhibition has two critical values: between these two critical PAGE 40 32 Value of Control Parameter Critical Value System Variable Value of Control Parameter Critical Value 2 Critical Value 1 System Variable values these systems are bistable, with two stable valu es separated by an unstable value. Changing the value of the parameter past either of these critical values will cause the variable to stabilize at one of the steady states, but this change in steady state is not permanent and can be reversed by changing the parameter value in the direction opposite to the original change (Fig. 13(b)). For each of these types of feedback the initial value of the variable will determine the steady state approached as the reaction occurs over time: small initials will tend towards the lower steady state while larger will tend towards the higher (Tyson, 2003). Figure 13. This figure shows the steady state curves that arise in positive feedback systems for a (a) one way switch and (b) two way switch. A simp le reaction network for the mutual inhibition type of regulation can be seen in Fig. 14: PAGE 41 33 S R A B Figure 14. This figure depicts the reaction network for mutual inhibition. asdfaaaaaaaaaa In this reaction network S promotes the production rate of R R th en increases the rate of conversion of A into B while A in turn decreases the amount of R by promoting its degradation rate, providing a positive feedback loop through a mutual inhibition mechanism. This system can be described with the following system of differential equations: d R [ ] dt = k 1 S k 2 R [ ] A [ ] k 2 R [ ] ( 2 31 ) d B [ ] dt = V m 1 A [ ] K m 1 + A [ ] R [ ] V m 2 B [ ] K m 2 + B [ ] ( 2 32 ) Under the assumption of conservation of mass: A [ ] = A 0 [ ] B [ ] ( 2 33 ) In these equations the concentrations of R and B are the dynamic variables under study and S is a c ontrollable parameter. Eq. ( 2 31 ) describes the dynamics of R The first term k 1 S ( ) in this equation describes the production of R as a zero order r eaction, and the last term k 2 R [ ] ( ) represents the constitutive degradation rate as a first order reaction. The middle term k 2 R [ ] A [ ] ( ) represents the degradation of R due to the feedback loop. The dynamics of B are giv en by Eq. ( 2 32 ) It is assumed that the conversions of A into B and PAGE 42 34 B into A are considered two different reactions, each governed by Michaelis Menten equations. When thi s system of equations is solved numerically the bistability can be easily observed with parameter values k 1 = 2, k 2 = 4, k 2 = 0.1 V m 1 = 2, V m 2 = 1 K m 1 = 1, and K m 2 = 1, and with S = 1.2, a value falling between the two critical values as seen in Fig. 13 (b). As seen in Fig. 15 the R concentration can converge to one of two steady states (both locally stable) corresponding to the upper or lower branches of the S shaped curve seen in Fig. 13(b). Where it converges depends strictly on the initi al concentrations of R and B Figure 15. Bistability arises in the mutual inhibition model. This plot shows the effect of initial concentrations of R around the middle branch of the S shaped curve seen in Fig. 13(b) in the numerical simulation. See te xt for selection of parameter values. Negative Feedback and Oscillation Another important regulatory mechanism is negative feedback, which can produce complex behaviors such as oscillations. There are two types of oscillating systems: damped and sustai ned. Damped oscillations can occur in negative feedback PAGE 43 35 S R A B loops with two or more components; these systems eventually settle on a stable steady state solution (Tyson, 2003). An example of this kind of network can be seen in Fig. 16. Figure 16. This figure shows the hypothetical reaction network that displays damping and sustained oscillations. In this reaction network, S promotes the degradation rate of R R increases the rate of conversion of A into B while A in turn increases the amount of R by promoting its production rate, completing a negative feedback loop. The equations describing a more complex system that can create oscillations are: d R [ ] dt = k 1 A [ ] k 2 S R [ ] k 2 R [ ] d B [ ] dt = V m 1 A [ ] K m 1 + A [ ] R [ ] V m 2 B [ ] K m 2 + B [ ] with the same conservation of mass equation from above (Eq. ( 2 33 ) ). Each component is the same as in the mutual inhibition, bistability example above. The same assumptions about the reactions governing R [ ] and the reactio ns for A [ ] and B [ ] made in the bistable system are made here. Damped oscillations can be seen in Fig. 17 under parameter values k 1 = 2, k 2 = 1, k 2 = 0.1 V m 1 = 6, V m 2 = 3 K m 1 = 0.5, K m 2 = 2, and S = 3 PAGE 44 36 Figure 17. Damped oscillations can occur in a negative feedback system. While initially the concentration of R alternates between higher and lower than the steady state, it eventually settles. See text for parameter values. Sustained oscilla tions about the steady state require a negative feedback loop with three or more components and proper parameters, where a mid reaction component causes a time delay in the feedback. Under the right parameter values the model used to create damping oscill ations can create the steady oscillations seen in Fig. 18. This behavior occurs when there are two or more possible steady states, and the concentration value will travel back and forth between them in the same pattern repeatedly, without deviation, until the reaction is either stopped or knocked out of this pattern via another method (Tyson, 2003). PAGE 45 37 Figure 18. Sustained oscillations can occur in a negative feedback system. The concentration of R alternates between two constants higher and lower than th e steady state, or between two steady states. Oscillations with more than two points can occur as well. The most complicated behavior possible is that of chaotic behavior. After the primary parameter of the model reaches a certain value in a nonlinear dynamical system with at least three dimensions, bifurcation and period doubling may occur; that is, when a time series simulation is run, the values that the product concentration of the reaction oscillates between will double, but a pattern will still be followed. As the value of that parameter increases, the amount of bifurcations will increase. However, after a certain parameter value is reached the reaction will no longer oscillate in a pattern between a fixed number of values; rather, it will start to jump between values in a pattern less, random behavior known as chaos. This state is very sensitive to the initial conditions of the system, such that even the smallest perturbation can have a huge effect on the output, and two initial values very clos e to each other can quickly evolve into two very different behaviors. This type of behavior can occur in either iterative or time series systems PAGE 46 38 involving either type of feed loop (feedback or feed forward), as long as the system is in some way nonlinear. PAGE 47 39 3. Methods 3.1. Model Discrimination Mathematical models have become one of the mo st important instruments in systems biology. However, multiple mathematical models describing the same reaction network with different interaction mechanisms may fit a set of experimental data equally well, and in order to determine what is actually happe ning in the network an ideal experiment must be designed to discriminate between the models by showing that only one model, but not the others, can fit optimized experimental data. The general form of the differential equation model for a reaction network is: x = dx dt = f x ( ) + g x ( ) u y = h x ( ) ( 3 1 ) where x is the state variable (concentrations) as an n dimensional vector, t is time, u is the input as a q dimensional vector, y is the output variable as an l dimensional vector, and f g and h are functions of the state variable, with g being a matrix sized n q In this for mulization it is assumed that the input u has no direct effect on the output y When a mathematical model is linear it can be written in the form x = A x + B u ( 3 2 ) y = Cx ( 3 3 ) where the matrices A B and C are constant matrices of sizes n n n q and l n respectively. When there is no change in the value of the state var iables over time (that is, they are constant), this is known as the steady state of the model, the value of which is PAGE 48 40 represented by x In Eqs. ( 3 1 ) and ( 3 2 ) this can occur when the time derivative is equal to zero: x = dx dt = 0 The steady state is also known as the equilibrium state or resting state For nonlinear models when u =0, such as that given in Eq. ( 3 1 ) the steady state can be calculated from solving the nonlinear system of equations: f x ( ) = 0 ( 3 4 ) Similarly for linear models, when u =0 such as Eq. ( 3 2 ) a steady state x satisfies A x = 0 If A is an invertible matrix then x = 0 In general, one way to study a dynamic system is to fo rce it out of its resting state via some sort of perturbation, and to study the behavior of its response. With the initial condition x 0 ( ) = x 0 the solution of Eq. ( 3 1 ) starting from x 0 has the general form x = x t x 0 ( ) ( 3 5 ) If this model has a steady state x then once the state vector in Eq. ( 3 5 ) is equal to this steady state value at time t 1 it remains equal to this point for all times t > t 1 A steady state is defined as locally asymptotically stable if the solutions represented by Eq. ( 3 5 ) converge to this value from every initial starting point close enough to the steady state. It is defined as globally asymptotically stable if the solutions converge to this value from every initial starting point. In other words x t x 0 ( ) x a s t for all x 0 near the steady state for local stability, and for all x 0 for global stability. PAGE 49 41 3.2. Symmetric Matrices and Quadratic Forms A real symmetric matrix W n n is a matrix whose transpose is equal to the original matrix: W T = W Here T represents the transpose of the matrix. Theorem : Any two eigenvectors of a symmetric matrix W n n are orthogonal. Proof : Let v 1 and v 2 be eigenvectors of W that correspond to distinct eigenvalues 1 and 2 It must be shown that v 1 T v 2 = 0 : 1 v 1 T v 2 = 1 v 1 ( ) T v 2 = W v 1 ( ) T v 2 = v 1 T W T ( ) v 2 = v 1 T W v 2 ( ) = v 1 T 2 v 2 ( ) = 2 v 1 T v 2 Hence 1 2 ( ) v 1 T v 2 = 0 and since 1 2 ( ) # 0 v 1 T v 2 = 0 and v 1 and v 2 are orthogonal. Theorem : A matrix W n n is orthogonally diagonalizable if and only if it is symmetric. Proof : Let A be an orthogonal matrix, where A 1 = A T and D diag onal matrix, such that W = A D A T ( 3 6 ) Taking the transpose of bo th sides of Eq. ( 3 6 ) W T = A D A T ( ) T = A T T D T A T = A D A T = W and thus W is symmetric (Lay, 2005). The second part of this proof, which is involved, will not be given here and can be found in (Horn and Johnson, 1985). Theorem : All eigenvalues of an n n symmetric matrix W n n are real. PAGE 50 42 Proof : Let denote any eigenvalue of W and v n the cor responding eigenvector such that W v = v ( 3 7 ) Define v n as the complex conjugate of v Then multiply both sides of Eq. ( 3 7 ) v T W v = v T v = v T v ( ) Here v T v is real a nd it must now be shown that v T W v is real (Sontag, 1998). Let q = v T W v where q To show that this is a real value it must be shown that q = q q q = v T W v = v T W v = v T W v = v T W v ( ) T = v T W T v = q Since both v T W v and v T v are real, must be real (Lay, 2005). For an n n symmetric matrix W n n a qua dratic form is a function V defined as V x ( ) = x T W x ( 3 8 ) where x is an n 1 vector (Dullerud, 2010). A quadratic form V is positive definite if V x ( ) > 0 for all x 0 and positive semidefinite if V x ( ) 0 for all x 0 These definitions carry over to the matrix W of the quadratic form: a positive definite matrix W is a symmetric matrix for which the quadratic form x T W x is positive definite, and similarly for positive semidefi nite. Theorem : For a symmetric matrix W n n let M be defined as M = m a x x T W x : x = 1 { } ( 3 9 ) The quadratic function given in Eq. ( 3 8 ) takes its maximum value as M at the eigenvector of W corresponding to the largest eigenvalue of W PAGE 51 43 Pro of : Since W is symmetric any two eigenvectors are orthogonal to each other and the matrix is orthogonally diagonalizable. L et M as defined in Eq. ( 3 9 ) be the greatest ei genvalue i of W Orthogonally diagonalize the matrix W as ADA T This can be achieved through the change of variables x = A y such that x T W x = y T D y ( 3 10 ) Also, x = A y = y for all y because A T A = I and A y 2 = y 2 and y = 1 if and only if x = 1 Here x T W x and y T D y have the same set of values as x and y range over the set of all unit vectors. D is the diagonal matrix with the eigenvalues of W a b c along the diagonal. To simplify the proof, assume W is a 3 3 matrix such that A has the eigenvectors u i as columns. The matrices D and A can be written as D = a 0 0 0 b 0 0 0 c # # # $ % & & & and A = [ u 1 u 2 u 3 ] If y is a vector with coordinates y 1 through y 3 then y T D y can be found to equal a : PAGE 52 44 y T D y = a y 1 2 + b y 2 2 + c y 3 2 a y 1 2 + a y 2 2 + a y 3 2 = a y 1 2 + y 2 2 + y 3 2 ( ) = a y 2 = a Thus M a by the definition of M However, y T D y equals a when y = e 1 = [ 1 0 0 ] T so M = a and from Eq. ( 3 10 ) the x that corresponds to e 1 is the eigenvector u 1 shown by x = A e 1 = u 1 u 2 u 3 [ ] 1 0 0 # # # $ % & & & = u 1 Thus M = e 1 T D e 1 = u 1 T W u 1 proving M and the quadratic form, is at its maximum when x is the eigenvector corresponding to the largest eigenvalue of W (Lay, 2005). The L 2 norm of the vector function y from a to b is defined as y 2 = y t ( ) T y t ( ) a b dt ( ) 1 2 ( 3 11 ) This will be used later. 3.3. Observability and the Observability Gramian In this section, two concept s from control theory will be summarized: observability and the observability gramian. They will be used later on during the derivation of the optimal experimental design. Observability is the ability to determine the initial state of a system from the o utput state. This concept assumes that the values of all constants and the state equation are known. A system model described by Eqs. ( 3 2 ) and ( 3 3 ) is defined as observable if the initial state x 0 ( ) can be determined from a known output y and a given input u over the time frame 0 t 1 [ ] where t 1 > 0 PAGE 53 45 Solving the linear differential Eq. ( 3 2 ) with x 0 = x 0 ( ) gives x t ( ) = e A t x 0 ( ) + e A t ( ) B u ( ) d 0 t # Inserting this into Eq. ( 3 3 ) gives y t ( ) = Ce A t x 0 ( ) + C e A t ( ) B u ( ) d 0 t # ( 3 12 ) Since y and u are known, this can be rewritten as Ce A t x 0 ( ) = y t ( ) ( 3 13 ) where y t ( ) = y t ( ) C e A t ( ) B u ( ) d 0 t # is known. Therefore, all that must be found is x 0 ( ) in Eq. ( 3 13 ) If u 0 then y t ( ) from Eq. ( 3 13 ) reduces to only y t ( ) = Ce A t x 0 ( ) The system is observable if and only if x 0 ( ) can be found from the measured output y t ( ) over a restricted time interval with zero input (Chen, 1984). The uniqueness of this initial condition x 0 ( ) depends on the singularity of the matrix Ce A t If it is nonsingular, that is, invertible, then the solution is unique; a singular matrix does not provide a unique solution. Above it is stated that the system model is observable if the initial state x 0 ( ) can be determined f rom a known output y and a given input u Theorem : A system model described by Eqs. ( 3 2 ) and ( 3 3 ) is observable if and only if there is a real, nonsingular n n matrix PAGE 54 46 W O t ( ) = e A T C T Ce A d 0 t ( 3 14 ) for all t > 0 Proof : The matrix function in Eq. ( 3 14 ) is found by pre multiplying from the left both sides of the Eq. ( 3 13 ) by e A T t C T and integrating over the time interval 0 t 1 [ ] to get e A T t C T Ce A t dt 0 t 1 # $ % & x 0 ( ) = e A T t C T y t ( ) dt 0 t 1 As lo ng as W O t 1 ( ) is nonsingular, a unique value for x 0 ( ) can be solved for: x 0 ( ) = e A T t C T Ce A t dt 0 t 1 # $ % & ( 1 e A T t C T y t ( ) dt 0 t 1 = W O ( 1 t 1 ( ) e A T t C T y t ( ) dt 0 t 1 This yields a unique solution to x 0 ( ) and therefore if W O t ( ) is no nsingular for any t > 0 the system model described by Eqs. ( 3 2 ) and ( 3 3 ) is observable. Now it must be shown that if W O t 1 ( ) is singular for all t 1 then the system model is not observable and no unique solution for x 0 ( ) can be found. If W O t 1 ( ) is singula r, there is an n 1 nonzero constant vector v that satisfies v T W O t 1 ( ) v = v T e A T C T Ce A v 0 t 1 d = Ce A v 2 2 d = 0 0 t 1 This implies Ce A t v 0 for all values of t in 0 t 1 [ ] From this, two initial states can be found that s olve a simplified version Eq. ( 3 12 ) when u =0: y t ( ) = Ce A t x i 0 ( ) 0 PAGE 55 47 These are x 1 0 ( ) = v 0 and x 2 0 ( ) = 0 and since there are two st ates that provide the same output, the solutions are not unique and the system is not observable. Hence, in order for the system model to be observable and for there to be a unique solution for the initial state the matrix W O t ( ) must be nonsingular, and the proof is complete. It should be noted that if a system is observable, both W O t ( ) and the matrix Ce A t are nonsingular and can be used to find a unique solution for the initial states for all values of t > 0 with a known output y and input u As can be seen from the derivations above, observability is dependent only on the constant matrices A and C and has no reliance on B Theorem : If all the eigenvalues of A have n egative real parts, that is, A is Hurwitz, then the unique solution of the equation A T W O + W O A = C T C ( 3 15 ) is the positive definite matrix W O = e A T C T Ce A d 0 # ( 3 16 ) known as the observability gramian. Proof : Insert Eq. ( 3 16 ) into Eq. ( 3 15 ) : A T W O + W O A = A T e A T C T Ce A + e A T C T Ce A A ( ) 0 # d = d dt 0 # e A T C T Ce A t ( ) d = e A T C T Ce A ( ) 0 = $ C T C This integral will converge only if A is Hurwitz. Thus the theorem is true. Theorem : The following two statements are equivalent: PAGE 56 48 1. The following matrix is of size n n and nonsingular for any t > 0 : W O t ( ) = e A T C T Ce A d 0 t 2. If all the eigenvalues of A have negative real parts, that is, A is Hurwitz, then the unique solution of the equation A T W O + W O A = C T C ( 3 17 ) is the positive definite matrix W O = e A T C T Ce A d 0 # ( 3 18 ) known as the observability gramian. Proof : If A is Hurwitz, then the unique solution of Eq. ( 3 17 ) can be expressed as in Eq. ( 3 18 ) The observability gramian is always positive semidefinite; it is positive definite if and only if it is nonsingula r. This establishes the equivalence of these two statements (Chen, 1999) When looking at the size of the output of the system, which is defined by the L 2 norm, after the system is released from x 0 = x 0 ( ) with u =0 and for a large enough t such that the steady state is reached, then y 2 2 = x 0 T e A T t C T Ce A t x 0 dt = 0 x 0 T e A T t C T Ce A t dt x 0 0 This can be simplified to: y 2 2 = x 0 T W O x 0 ( 3 19 ) Where W O is the observability gramian defined as previously. Here the sizes of the eigenvalues of the observability gramian describe the L 2 norm of the output of the model PAGE 57 49 (Unneland, 2007). Since Eq. ( 3 19 ) is in quadratic form, this form will take its maximum value at the eigenvector corresponding to the largest eigenvalue. 3.4. Multiple Models The purpose of the following sections will be to set up the process of differentiating between two models for the same network. It is assumed that each model has the same number of state variables, the same steady states, and fit the same data equally well. The first step is to link the two models, then to find the difference between their outputs. x 1 x 2 # $ % & = x = f 1 x 1 ( ) f 2 x 2 ( ) # # $ % & & = f x ( ) # $ % & % + g 1 x 1 ( ) g 2 x 2 ( ) # # $ % & & = g x ( ) # $ % & % u y = y 1 y 2 = h 1 x 1 ( ) h 2 x 2 ( ) = h x ( ) # $ % % & % % Here x i represents state variables of model i f i h i and g 1 are functions of x i y i is the output of model i and u is the shared input for both models ( i =1, 2 ). An op timal experiment can be defined as an experiment that maximizes the value of the difference between the outputs y = y 1 y 2 over a set of allowable experimental perturbations in the initial state conditions. The specific difference this me thod attempts to maximize is the L 2 distance between the outputs, or O pt i m al E x pe r i m e nt = m a x y 2 { } = m a x y 1 y 2 2 { } Assumptions are made that the steady states of the two models are asymptotically stable and the same. When the mathematical models are linear the two models can be written in the compact form PAGE 58 50 x = A x + B u ( 3 20 ) y = Cx ( 3 21 ) where x = x 1 x 2 # $ % & A = A 1 0 0 A 2 # $ % & B = B 1 B 2 # $ % & C = C 1 C 2 [ ] ( 3 22 ) where the matrices A B and C are redefined as the constant matrices of sizes 2 n 2 n 2 n q and l 2 n respectively. The matrices A 1 and A 2 are Hurwitz; that is, their eigenvalues have negative real parts so they guarantee that the steady states are asymptotically stable: Re i [ ] < 0 f or a l l i To study the behavior of two solutions during the return to the shared stable steady state, optimized initial conditions must be found such that x 1 0 ( ) = x 2 0 ( ) maximizes y 2 for the two models when u =0. The optimum initial condi tions can be found by first finding a positive semidefinite matrix P that will be shown to be equivalent to the observability gramian (Mlykti, 2010). To find P the definition of the L 2 distance norm is used. 3.5. Derivation of P from L 2 Norm For the sys tem Eqs. ( 3 20 ) and ( 3 21 ) the solution for x t ( ) with initial con dition x 0 ( ) = x 0 when u =0, is x t ( ) = e A t x 0 ( 3 23 ) Substituting Eq. ( 3 23 ) into the output y t ( ) in Eq. ( 3 21 ) gives: PAGE 59 51 y t ( ) = Cx t ( ) = Ce A t x 0 ( 3 24 ) The L 2 norm squared of y is defined as: y 2 2 = y t ( ) T y t ( ) 0 dt ( 3 25 ) Placing Eq. ( 3 24 ) into Eq. ( 3 25 ) considering x 0 is a constant vector, gives y 2 2 = x 0 T e A T t C T Ce A t x 0 dt = 0 x 0 T e A T t C T Ce A t dt x 0 0 ( 3 26 ) From Eq. ( 3 26 ) the matrix function P 2 n 2 n is d efined as P = e A T t C T Ce A t dt 0 ( 3 27 ) Then Eq. ( 3 26 ) simplifies to the quadratic form y 2 2 = x 0 T P x 0 ( 3 28 ) where the initial condition vector x 0 has the form of: y = y 1 y 2 a nd x 0 = x 1 T 0 ( ) x 2 T 0 ( ) # $ % T ( 3 29 ) Here P is equivalent to the observability gramian. As can be seen here, the difference between the two model outputs is determined by the observability matrix P 3.6. Optimum Initial Conditions An optimized experimental de sign would have a single set of initials for both models, therefore the initial condition in Eq. ( 3 29 ) must be in the form x 0 = x x # $ % & where x 1 0 ( ) = x 2 0 ( ) = x To make this condition true P is partitioned into blocks: PAGE 60 52 P = P 11 P 12 P 12 T P 22 # $ % & ( 3 30 ) where each block is of size n n In order to show that P can be separated in this fashion Eq. ( 3 24 ) must first be written as such: y t ( ) = C 1 C 2 [ ] e A 1 0 0 A 2 # $ % & t x 1 0 ( ) x 2 0 ( ) # $ % & ( 3 31 ) This equation uses the definitions of the matrices A and C g iven in Eq. ( 3 22 ) When this equation for y is put into Eq. ( 3 25 ) and expanded the fo llowing can be seen: y 2 2 = x 1 T e A 1 T t C 1 T C 1 e A 1 t x 1 x 1 T e A 1 T t C 1 T C 2 e A 2 t x 2 x 2 T e A 2 T C 2 T C 1 e A 1 t x 1 + x 2 T e A 2 T t C 2 T C 2 e A 2 t x 2 # $ % & 0 ( ) dt ( 3 32 ) When each of the four components of this equation are separated this simplifies to y 2 2 = x 1 T P 11 x 1 + x 1 T P 12 x 2 + x 2 T P 12 T x 1 + x 2 T P 22 x 2 ( 3 33 ) wher e P 11 = e A 1 T t C 1 T C 1 e A 1 t 0 dt P 12 = e A 1 T t C 1 T C 2 e A 2 t 0 # dt P 12 T = e A 2 T t C 2 T C 1 e A 1 t 0 # dt P 22 = e A 2 T t C 2 T C 2 e A 2 t 0 dt Note that P 12 T = P 21 Thus P can be partitioned as in Eq. ( 3 30 ) When x 1 0 ( ) = x 2 0 ( ) = x is true, Eq. ( 3 33 ) can be written is: y 2 2 = x 1 T P 11 x 1 + x 1 T P 12 x 2 + x 2 T P 12 T x 1 + x 2 T P 22 x 2 = x T P 11 + P 12 + P 12 T + P 22 ( ) x = x T R x ( 3 34 ) PAGE 61 53 where R = P 11 + P 22 + P 12 + P 12 T Here R is a real symmetric matrix and thus Eq. ( 3 34 ) is in quadratic form. It should be noted that the dimensions of P and x are of dimensions 2 n 2 n and 2 n 1 respectively, while R and x are of dimensions n n and n 1 respectively. Similar to Eq. ( 3 28 ) above, the value of Eq. ( 3 34 ) is the difference in the final outputs of each model in quadratic form, which serves to discriminate between each model and is largest when x is the eigenvect or corresponding to the maximum eigenvalue of the matrix R is used as the optimized initial state ( Mlykti, 2010 ). For non linear models this method can be applied by linearizing the models around a locally stable steady state and using the resulting Jacobian matrix as the matrix A 3.7. Application of Method to Enzyme Inhibition Models To show how this method works for any two models it has been applied to the well known case of competitive versus uncompetitive inhibition in en zyme kinetics. First the two rivaling mechanisms will be looked at, and their mathematical models will be derived. Competitive Inhibition The competitive inhibition mechanism shown below and the associated model will be referred to as mechanism one, with unique constants k i 1 k ni 1 and k n 21 PAGE 62 54 k 1 k 2 E + S + I ES E + P EI k n 1 k n 21 k ni1 k i1 Here the variables E S ES EI and P are defined as in Chapter 2. The total concentration equations for the enzyme, substrate, and inhibitor are: E 0 # $ = E 1 # $ + E S 1 # $ + E I 1 # $ I 0 # $ = I 1 # $ + E I 1 # $ S 0 # $ = S 1 # $ + E S 1 # $ + P 1 # $ The three dimensional differential equation model of this mechanism can be written as d E 1 # $ dt = % k 1 E 1 # $ S 1 # $ + k n 1 + k 2 ( ) E S 1 # $ % k ni 1 E 1 # $ I 1 # $ + k i 1 E I 1 # $ % k n 21 E 1 # $ P 1 # $ d S 1 # $ dt = % k 1 E 1 # $ S 1 # $ + k n 1 E S 1 # $ d I 1 # $ dt = % k ni 1 E 1 # $ I 1 # $ + k i 1 E I 1 # $ ( 3 35 ) From the conser vation equations E I 1 [ ] E S 1 [ ] and P 1 [ ] can be written as E I 1 # $ = I 0 # $ % I 1 # $ E S 1 # $ = E 0 # $ % E 1 # $ % I 0 # $ + I 1 # $ P 1 # $ = S 0 # $ % S 1 # $ % E 0 # $ + E 1 # $ + I 0 # $ % I 1 # $ Uncompetitive Inhibition The uncompetitive inhibition mechanism shown below and the associated model wi ll be referred to as mechanism two, with unique constants k i 2 k ni 2 and k n 22 PAGE 63 55 E + S ES + I E + P k 1 k n 1 k 2 ESI k n22 k i 2 k ni2 Here the variable ESI is defined as in Chapter 2. The total concentration equations for the enzyme, substrate, and inhibitor are: E 0 # $ = E 2 # $ + E S 2 # $ + E SI 2 # $ I 0 # $ = I 2 # $ + E SI 2 # $ S 0 # $ = S 2 # $ + E S 2 # $ + E SI 2 # $ + P 2 # $ The three dim ensional differential equation model of mechanism two can be written as d E 2 # $ dt = % k 1 E 2 # $ S 2 # $ + k n 1 + k 2 ( ) E S 2 # $ % k n 22 E 2 # $ P 2 # $ d S 2 # $ dt = % k 1 E 2 # $ S 2 # $ + k n 1 E S 2 # $ d I 2 # $ dt = % k ni 2 E S 2 # $ I 2 # $ + k i 2 E SI 2 # $ ( 3 36 ) From the conservation equations E SI 2 [ ] E S 2 [ ] and P 2 [ ] can be written as E SI 2 # $ = I 0 # $ % I 2 # $ E S 2 # $ = E 0 # $ % E 2 # $ % I 0 # $ + I 2 # $ P 2 # $ = S 0 # $ % S 2 # $ % E 0 # $ + E 2 # $ Steady State Calculations Now the steady states for all three components, subst rate, enzyme, and inhibitor, of each model must be equal. The steady states for the enzyme, substrate, and inhibitor PAGE 64 56 of model 1 are defined as E 1 # $ S 1 # $ and I 1 # $ respectively, and can be solved for by setting the differential equations equal to zero: d E 1 [ ] dt = 0 d S 1 [ ] dt = 0 a nd d I 1 [ ] dt = 0 The steady states of model 2 are found similarly and defined as E 2 # $ S 2 # $ and I 2 # $ Now the parameters must be selected such that: E 1 # $ = E 2 # $ S 1 # $ = S 2 # $ I 1 # $ = I 2 # $ Linearization The linearized system of equations corresponding to the six dimensional nonlinear system from Models 1 and 2, from Eqs. ( 3 35 ) and ( 3 36 ) takes the form x A x = as: E 1 S 1 I 1 E 2 S 2 I 2 # # # # # # # # $ % & & & & & & & & = x 1 x 2 x 3 0 0 0 x 4 k 1 E 1 $ % k n 1 0 0 0 k ni 1 I 1 $ % 0 x 5 0 0 0 0 0 0 y 1 y 2 y 3 0 0 0 y 4 k 1 E 2 $ % k n 1 0 0 0 k ni 2 I 2 $ % 0 y 5 # # # # # # # # # # $ % & & & & & & & & & & E 1 S 1 I 1 E 2 S 2 I 2 # # # # # # # # $ % & & & & & & & & where PAGE 65 57 x 1 = k n 1 k 2 k n 21 S 0 # $ % + k n 21 E 0 # $ % k n 21 I 0 # $ % 2 k n 21 E 1 # $ % + k n 21 k ni 1 ( ) I 1 # $ % + k n 21 k 1 ( ) S 1 # $ % x 2 = k n 21 k 1 ( ) E 1 # $ % x 3 = k n 21 k ni 1 ( ) E 1 # $ % + k n 1 + k 2 k i 1 x 4 = k 1 S 1 # $ % k n 1 x 5 = k ni 1 E 1 # $ % k i 1 y 1 = k n 1 k 2 k n 22 S 0 # $ % + k n 22 E 0 # $ % + k n 22 k 1 ( ) S 2 # $ % 2 k n 22 E 2 # $ % y 2 = k n 22 k 1 ( ) E 2 # $ % y 3 = k n 1 + k 2 y 4 = k 1 S 2 # $ % k n 1 y 5 = k ni 2 E 2 # $ % 2 k ni 2 I 2 # $ % + k ni 2 I 0 # $ % k ni 2 E 0 # $ % k i 2 The paramet ers used to compare the models are fixed at the values shown in Table 2. Table 2. The values of parameters for optimization of initial conditions for inhibition models. Parameter k 1 k n 1 k 2 k n 21 k n 22 k i 1 k ni 1 k i 2 k ni 2 E 0 [ ] S 0 [ ] I 0 [ ] Value 10 500 0.1 0.001 0. 001006 10 1 6.6498 1 1 100 10 From above the steady states are calculated as E 1 # $ = E 2 # $ = 0.3805 S 1 # $ = S 2 # $ = 33.2490 I 1 # $ = I 2 # $ = 9.6335 The 6 6 matrix above is used as the matrix A in determining P and C is the unit matrix. It sh ould be noted that C can be altered so that only one variable may be studied: by setting the columns corresponding to some of the variables as a zero vector, the analysis of a chosen variable, corresponding to an unaltered vector, can be achieved. Howeve r, here C will be the full unit matrix to study all the chosen variables: PAGE 66 58 C = C 1 C 2 [ ] = 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 # $ $ $ % & The matrix P is then partitioned, the matrix R is found, and the eigenvector corresponding to the largest eigenvalue is used in conjunction with the calcu lated values of the steady states as the optimal initial states to maximize the difference in the models' outputs. The optimal initial conditions are calculated as E I ni t [ ] = 0.3823 S I ni t [ ] = 33.2472 I I ni t [ ] = 8.6334 using the method detailed above. The figures below demonstrate the effectiveness of this method. Fig. 19 shows the difference between the models as a function of time for two sets of arbitrary initials in (a) and the optimized initials in (b). As can be seen, when starting the system at any non optimal initial condi tion the difference between the models will always be smaller than when the system is started at optimized initials. (a) PAGE 67 59 (b) Figure 19. Difference between the outputs of Competitive and Uncompetitive Inhibition Models. This figure displays the res ults for the difference between the models over time when started from (a) two sets of arbitrary initials around the steady state and (b) optimized initials E I ni t [ ] = 0.3823 S I ni t [ ] = 33.2472 and I I ni t [ ] = 8.6334 PAGE 68 60 C D A B S C D A B S 4. Results Starting with only experimental data, this section will demonstrate how to use the method to discriminate between two positive feedback mechanisms of a hypothetical reaction network. Then a five step experimental protocol for discrimination between rivaling mechanisms for any given reaction network when starting from only experimental data will be detailed. 4.1. Differentiation o f Two Feedback Mechanisms from Experimental Data The competing feedback mechanisms are shown in Fig. 20. (a) (b) Figure 20. This figure shows the two distinct positive feedback mechanisms. (a) shows mechanism 1, where the diamond indicates that the feedback reduces the rate of conversion. (b) shows mechanism 2; the large arrow indicates an increase in conversion rate. In this network letters A B C and D represent the components of the network. In mechanism 1 seen in Fig. 20(a), the feedback reduces the conversion rate of B back into A ; in mechanism 2, seen in Fig. 20(b), the feedback increases the conversion rate of A into B The following set of two differential equations are used to model mechanism 1: PAGE 69 61 d B [ ] dt = k 1 S [ ] A [ ] k 2 B [ ] 1 + 1 D [ ] d D [ ] dt = k 3 C [ ] B [ ] k 4 D [ ] ( 4 1 ) The first equation describes the dynamics of B [ ] The first term of this equation describes the conversion rate of A into B as a first order reaction with rate constant k 1 S [ ] The second term represents the rate of conversion of B back into A The denominator of the se cond term 1 + 1 D [ ] ( ) of the first equation describes the effect of the feedback in this mechanism, where 1 is the parameter that measures the strength of the feedback. Here the feedback decreases the conversion r ate of B back into A which is k 2 B [ ] in the absence of feedback. [S] is defined as a parameter required for the reaction to begin. The second equation describes the dynamics of D [ ] where the conversion rate of C into D as a second order reaction with rate constant k 3 and the reverse reaction of D back into C as a first order reaction with rate constant k 4 The following set of two differential equations are us ed to model mechanism 2: d B [ ] dt = S [ ] A [ ] k 1 + 2 D [ ] ( ) k 2 B [ ] d D [ ] dt = k 3 C [ ] B [ ] k 4 D [ ] ( 4 2 ) The first equation again desc ribes the dynamics of B [ ] The first term describes the conversion rate of A into B Here the positive feedback augments the forward rate constant k 1 for A into B by 2 D [ ] where 2 i s the strength of the feedback. The second term represents the conversion of B back into A as a first order reaction with rate constant k 2 The second equation is the same as the second equation in Eq. ( 4 1 ) Note that when PAGE 70 62 1 = 2 = 0 both models become the same. Both models assume the following equations are true: A [ ] = A 0 [ ] B [ ] C [ ] = C 0 [ ] D [ ] where A 0 [ ] and C 0 [ ] are the initial concentrations of A and C respectively. The steady states B 1 # $ and D 1 # $ of model 1 can be calculated by setting the time derivatives in Eq. ( 4 1 ) equal to zero, namely d B [ ] dt = d D [ ] dt = 0 : B 1 # $ = % 1 + k 1 k 3 A 0 [ ] S [ ] + & 1 2 k 3 k 1 S [ ] + k 2 + k 1 1 C 0 [ ] S [ ] ( ) D 1 # $ = % 1 ( k 1 k 3 A 0 [ ] S [ ] + & 1 2 k 1 1 S [ ] k 3 A 0 [ ] + k 4 ( ) where 1 = k 1 k 3 1 A 0 [ ] C 0 [ ] S [ ] # k 2 k 4 # k 1 k 4 S [ ] and 1 = k 1 2 k 3 2 A 0 [ ] 2 S [ ] 2 + 2 k 1 2 k 3 2 1 A 0 [ ] 2 C 0 [ ] S [ ] 2 + 2 k 1 k 2 k 3 k 4 A 0 [ ] S [ ] + 2 k 1 2 k 3 k 4 A 0 [ ] S [ ] 2 + k 1 2 k 3 2 1 2 A 0 [ ] 2 C 0 [ ] 2 S [ ] 2 # 2 k 1 k 2 k 3 k 4 1 A 0 [ ] C 0 [ ] S [ ] + 2 k 1 2 k 3 k 4 1 A 0 [ ] C 0 [ ] S [ ] 2 + k 2 2 k 4 2 + 2 k 1 k 2 k 4 2 S [ ] + k 1 2 k 4 2 S [ ] 2 Similarly for model 2 the steady states B 2 # $ and D 2 # $ are: B 2 # $ = % 2 + k 1 k 3 A 0 [ ] S [ ] + & 2 2 k 3 k 1 S [ ] + k 2 + 2 C 0 [ ] S [ ] ( ) D 2 # $ = % 2 ( k 1 k 3 A 0 [ ] S [ ] + & 2 2 2 S [ ] k 3 A 0 [ ] + k 4 ( ) PAGE 71 63 where 2 = k 3 2 A 0 [ ] C 0 [ ] S [ ] # k 2 k 4 # k 1 k 4 S [ ] and 2 = k 3 2 2 2 A 0 [ ] 2 C 0 [ ] 2 S [ ] 2 + 2 k 1 k 3 2 2 A 0 [ ] 2 C 0 [ ] S [ ] 2 + 2 k 1 k 3 k 4 2 A 0 [ ] C 0 [ ] S [ ] 2 # 2 k 2 k 3 k 4 2 A 0 [ ] C 0 [ ] S [ ] + k 1 2 k 3 2 A 0 [ ] 2 S [ ] 2 + 2 k 1 2 k 3 k 4 A 0 [ ] S [ ] 2 + 2 k 1 k 2 k 3 k 4 A 0 [ ] S [ ] + k 1 2 k 4 2 S [ ] 2 + 2 k 1 k 2 k 4 2 S [ ] + k 2 2 k 4 2 Notice that when 2 = k 1 1 ( 4 3 ) then the steady states of each model become equal: B 1 # $ = B 2 # $ and D 1 # $ = D 2 # $ ( 4 4 ) 4.2. Generation of Artificial Experimental Data Under the condition in Eq. ( 4 3 ) both models have the same steady state for any set of parameter values. The parameter values are now fixed at the values in Table 3. Table 3. The parameter values for the two positive feedback models for opti mal initial condition design. Parameter k 1 k 2 k 3 k 4 S [ ] 1 2 A 0 [ ] C 0 [ ] Value .01 5 3 1.5 35 .12 .0012 10 10 The steady states for this set of parameter values are calculated as: B 1 # $ = B 2 # $ = 1.1436 D 1 # $ = D 2 # $ = 6.9580 ( 4 5 ) Mechanism 2 is assumed to be the correct mechanism, and since no actual data is available for this study, model 2 is used to generate artificial experimental data. The model is solved numerically with the initial conditions B 0 [ ] = 0 and D 0 [ ] = 0 PAGE 72 64 and time series data for B [ ] and D [ ] is generated. An error with a 5% standard deviat ion around the concentration values at each time point from a normal distribution is added to mimic experimental error. The artificially created data for B [ ] and D [ ] at each time point is listed in Table 4. Table 4. Data points artificially created for comparison of Models 1 and 2 from initial values B 0 [ ] = 0 and D 0 [ ] = 0 Time [B]+Error [D]+Error 0.0 0.0000 0.0000 0.2 0.4210 1.3340 0.4 0.6979 3.3404 0.6 0.9067 4.72 34 0.8 0.9724 5.7943 1.0 0.9378 6.6427 1.2 1.0529 6.3360 1.4 1.1563 6.3120 1.6 1.1235 6.3259 1.8 1.1636 7.0311 2.0 1.1719 6.8564 2.2 1.1490 7.0278 2.4 1.2455 6.9740 2.6 1.1207 6.5088 2.8 1.1004 7.2264 3.0 1.1590 7.0795 3.2 1.1548 7.0766 3.4 1 .1211 7.0968 3.6 1.1539 7.4366 3.8 1.1108 6.4247 4.0 1.1453 7.3518 4.2 1.1035 6.5412 4.4 1.1461 7.3105 4.6 1.1456 7.2041 4.8 1.0487 7.1154 5.0 1.1931 7.0367 Fig. 19 shows the simulation of the two models with the parameter values in Table 3 with initial conditions B 0 [ ] = 0 and D 0 [ ] = 0 The solid line represents model 1, PAGE 73 65 the dashed line is model 2, and the dots are the data given in Table 4 above. It can be seen that in terms of how well they fit the data, the models are virtually indistinguishable. PAGE 74 66 Figure 21. This figure shows the numerical solution of Models 1 and 2 compared with the artificially generated data. The initial conditions are B 0 [ ] = 0 and D 0 [ ] = 0 an d the values for the parameters are listed in Table 3. As seen, both models share a steady state and capture the data equally well. PAGE 75 67 4.3. Optimal Experimental Design Now an optimal initial condition, which maximizes the difference between the outputs of the two models, can be calculated using the method described in Chapter 3. The differential equations from the two mechanisms are now taken as one larger system: d B 1 [ ] dt = k 1 S [ ] A 0 [ ] B 1 [ ] ( ) k 2 B 1 [ ] 1 + 1 D 1 [ ] d D 1 [ ] dt = k 3 C 0 [ ] D 1 [ ] ( ) B 1 [ ] k 4 D 1 [ ] d B 2 [ ] dt = S [ ] A 0 [ ] B 2 [ ] ( ) k 1 + 2 D 2 [ ] ( ) k 2 B 2 [ ] d D 2 [ ] dt = k 3 C 0 [ ] D 2 [ ] ( ) B 2 [ ] k 4 D 2 [ ] where B 1 [ ] and D 1 [ ] are the variables in model 1, and similarly for B 2 [ ] and D 2 [ ] The linearized form of this system around the shared steady states in Eq. ( 4 4 ) is written in matrix form: B 1 D 1 B 2 D 2 # # # # # $ % & & & & & = x 1 x 2 0 0 x 3 x 4 0 0 0 0 y 1 y 2 0 0 y 3 y 4 # # # # $ % & & & & B 1 D 1 B 2 D 2 # # # # $ % & & & & where PAGE 76 68 x 1 = k 1 S [ ] k 2 1 + 1 D 1 # $ % & x 2 = k 2 1 B 1 # $ % & 1 + 1 D 1 # $ % & ( ) 2 x 3 = k 3 C 0 [ ] D 1 # $ % & ( ) x 4 = k 3 B 1 # $ % & k 4 y 1 = k 1 S [ ] 2 D 2 # $ % & S [ ] k 2 y 2 = 2 A 0 [ ] S [ ] 3 B 2 # $ % & S [ ] y 3 = k 3 C 0 [ ] D 2 # $ % & ( ) y 4 = k 3 B 2 # $ % & k 4 The method of initial condition design described in Chapter 3 is applied with the parameter values in Table 3, and the optimized initial conditions are calculated: B 0 [ ] = 0.1437 D 0 [ ] = 6.9685 A new set of data is artificially created from same model (model 2) with these initial conditions using the same method as above with the same parameter values Fig. 20 shows the simulation of the two models wi th the same parameter values in Table 3 starting from the optimal initial conditions. In this plot the solid line represents model 1, the dashed line is model 2, and the dots are the data from the new initial conditions. A comparison of Figs. 21 and 22 shows that the optimal initial conditions maximize the difference between the models' outputs as the outputs approach the shared steady states. When both models start at the zero initial conditions, they approach the steady states with the same behavior throughout the simulation. With the optimal initial conditions, the models still approach the steady states with similar behaviors, but model 2 PAGE 77 69 now reaches the steady states faster, and it is this difference that allows for the determination of the correc t model. Once the simulations have reached the steady states, the models cannot be discriminated as they both match the data equally. The data for both Figs. 21 and 22 was generated from model 2, with an error of 5%, to simulate experimental error. If e rrors of less than 5% are used to generate data then the data points would lie closer to the curves representing model 2, and with 0% error the data would exactly match the curve for model 2. PAGE 78 70 Figure 22. This figure shows a comparison of the time series of Models 1 and 2 with the generated artificial data when started from the optimum initials B 0 [ ] = 0.1437 and D 0 [ ] = 6.9685 As can be seen, model 2 captures the data better than model 1. PAGE 79 71 The difference between the variables of the two competing models and the effect of the optimal initial conditions can be seen more clearly if this difference is plotted over time. Fig. 21 shows this difference over time for (a) three runs with arbitrary initials and (b) the optimum initials B 0 [ ] = 0.1437 and D 0 [ ] = 6.9685 (b). As can be seen, the difference between the models when started from the optimum initials is always greater than the difference when started from arbitrary initials. The scale of the y axis has been made equal in both figures to make comparison easier. (a) PAGE 80 72 (b) Figure 21. This figure depicts the difference between the model outputs when started from (a) three sets of arbitrary initials and (b) optimized initials B 0 [ ] = 0.1437 and D 0 [ ] = 6.9685 As can be seen, the difference between the optimized initials is greater than that between any of the arbitrary initials. 4.4. Complete Experimental Protocol Based on the steps above, a five step experimental protocol to discriminate between to rival mechanisms can be summarized as: 1. Collect experimental time series data from a set of initial conditions. 2. Propose two competing mechanisms and estimate the parameters from this data set by minimizing the diff erence between the model output and the experimental data for each model. Ensure the steady states of the two models are equal. PAGE 81 73 3. Calculate the optimal initial conditions that maximize the difference between the model outputs. 4. Run a second experiment (optimized) to generate a second time series data set. 5. Simulate the two competing models starting from the optimized initials found in step 3 and compare the model results with the second experimental data set. PAGE 82 74 5. Discussion and Conclusion In this study a method for designing an optimal experiment that allows for the discrimination of multiple competing mechanisms for a biological network has been described. The specific method studied was optimal initial condition design for the state variables, and it was demonstrated on two different reaction networks and allowed for the invalidation of incorrect models This study began with a description of the background of systems biology and a problem that occurs when attempting to model the behavior of a biologic al system; namely, two competing reaction mechanisms that capture experimental data equally Next the details of mathematically modeling reaction networks were given, and both simple and complex behaviors were studied. Mathematical methods of studying re action rates and steady states were described, and the biology of enzyme inhibition kinetics was discussed. Then the necessary theory needed for optimal initial condition design was given, followed by the method itself. The use of the method on competing versions of the enzyme inhibition mechanisms, specifically competitive and uncompetitive inhibition, was presented as an example of the application of the method. Finally, a network with two competing, hypothetical, positive feedback mechanisms was used to demonstrate the complete experimental protocol starting from only experimental data. The method was first used on an example concerning enzyme inhibition kinetics as proof that the method works. This example was chosen to demonstrate the effectiveness of the method due to its well known behaviors that standard methods use for discrimination. With traditional discrimination methods, examination of the PAGE 83 75 Michaelis Menten equations and Lineweaver Burk plots would be used to differentiate these models by cha nging the inhibitor and initial substrate concentrations and focusing on the behavior of the reaction as it progresses. These methods operate under certain assumptions, such as the equilibrium assumptions on inhibitor binding, and would require multiple r uns of the experiment with varying concentrations of the substrate and inhibitor, which could be both time consuming and expensive. Better methods to discriminate between these mechanisms would require fewer assumptions and fewer runs to determine the cor rect mechanism. In this study the only assumptions made about the mechanisms were that the steady states and initial conditions of the variables for each model were the same. The parameters were chosen such that the steady states were forced to be equal a nd that from arbitrary initials the simulations of each mechanism were nearly identical. Through the method of optimal initial condition design it was shown that the differences between the model outputs could be maximized such that discrimination could b e achieved with only one more experimental run; therefore the method works as described. The second network consisted of two theoretical mechanisms with positive feedback loops: one that would increase the forward conversion rate of the first reaction, an d one that would decrease the reverse conversion rate of the same reaction. In this example the method was used when starting from only experimental data to detail the experimental protocol. As in the enzyme inhibition example, parameters were chosen suc h that the mechanisms shared a steady state and time series from many arbitrary initials would fit the data equally well. It was found that the method of optimization increased the differences between the models during the time frame as the steady state PAGE 84 76 w as approached: the models still showed similar behavior over the course of the simulation, but one model would reach the final value much later than the other model. The feedback mechanisms were also used to give a step by step account of the entire exper imental protocol. Once the initial experiment has been run, the parameters of the models are estimated from the data once the steady states of each competing mechanism are forced to be equal. Once every parameter is given a value and the shared steady st ates are found, the optimal initial conditions are calculated. The experiment is run again with the newly found initials and data is plotted. The two models are simulated with the new initials, and via comparison with the data, one model is invalidated. In this study, the data used to demonstrate the protocol was artificially created from a model and known set of parameters. Under the standard experimental protocol the data would be available at the start of the process, with the models and parameters fi t to the data such that the steady states are equal. However, for this study no such data was readily available, and the results from such data could be easily replicated via the artificial data. Areas of future work include expanding the method to includ e comparing more than two mechanisms or studying a greater number of variables at once, and to design optimal experiments based on other components of the model, such as the input or parameters. The method detailed in this study was designed such two comp eting mechanisms could be discriminated between, and was applied to enzyme inhibition and feedback examples with two and three state variables, respectively. However, this method could be applied to networks with any number of variables, and, through mani pulation of the C matrix of the model output equation, specific variables can be PAGE 85 77 selected for or removed from analysis. It would also be beneficial to be able to discriminate between more than two competing mechanisms at a time. While this could be achie ved by comparing two models at a time for every involved mechanism, it would be more helpful to alter the method such that more than two mechanisms can be compared simultaneously. Methods to discriminate between models, and thus corresponding experimental protocols, could be derived based on optimal input design, or on structural changes based on numerical changes of the parameters. PAGE 86 78 Bibliography [1] Berg, J., Tymoczko, J., & Stryer, L. (2006). Biochemistry (Sixth Edit., p. 1026). New York: W. H. Free man and Company. [2] Bisswanger, Hans (2008). Enzyme Kinetics Principles and Methods (Second Revised and Updated Edit.). Weinheim, Germany: Wiley VCH Verlag GmbH & Co. KGaA. [3] Chen, C. T. (1984). Linear System Theory and Design New York: CBS C ollege Publishing. [4] Chen, C. T. (1999). Linear System Theory and Design Electrical Engineering (Third Edit., p. 348). New York: Oxford University Press, Inc. [5] Csiksz Nagy, A., Kapuy, O., Tth, A., Pl, C., Jensen, L. J., Uhlmann, F., et al. (2009). Cell cycle regulation by feed forward loops coupling transcription and phosphorylation. Molecular systems biology 5 6. [6] Dullerud, G. E., & Paganini, F. G. (2010). A Course in Robust Control Theory: A Convex Approach International Jour nal of Robust and Nonlinear Control Springer Science+Business Media, Inc. [7] Horn, R. A. & Johnson, C. R. (1985). Matrix Analysis Cambridge, New York: Cambridge University Press. [8] Lay, David C. (2005). Linear Algebra and Its Applications (Th ird Updated Edit.). Reading, Massachusetts: Addison Wesley. [9] Leskovac, Vladimir (2003). Comprehensive Enzyme Kinetics New York: Kluwer Academic/Plenum Publishers. [10] Mlykti, B., August, E., Papachristodoulou, A., & El Samad, H. (2010). Di scriminating between rival biochemical network models: three approaches to optimal experiment design. BMC systems biology 4 (38), 16. [11] Mitrophanov, A. Y., & Groisman, E. A. (2008). Positive feedback in cellular control systems. BioEssays: news a nd reviews in molecular, cellular and developmental biology 30 (6), 542 55. PAGE 87 79 [12] Sontag, E. D. (1998). Mathematical Control Theory: Deterministic Finite Dimensional Systems (Second Edi., p. 531+xvi). New York: Springer Verlag. [13] Tyson, J., Chen, K. C., & Novak, B. (2003). Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Current Opinion in Cell Biology 15 (2), 221 231. [14] Unneland, K., Van Dooren, P., & Egeland, O. (2007). A Novel Scheme for P ositive Real Balanced Truncation. 2007 American Control Conference 947 952. Programs Matlab The Language of Technical Computing: www.mathworks.com Maple The Essential Tool for Mathematics and Modeling: www.maplesoft.com 