From Sequence to Function Through Secondary Structure Kinetics of RNA and DNA Citation Meruelo, Alejandro Daniel (2006) From Sequence to Function Through Secondary Structure Kinetics of RNA and DNA. Senior thesis (Major), California Institute of Technology. doi:10.7907/4bxt-wa89. https://resolver.caltech.edu/CaltechTHESIS:02132018-142852598 Abstract A number of efforts to determine function from sequence of RNA and DNA have been made with varying success. Here we study the determination of function from sequence of DNA and RNA through their secondary structure kinetics, specifically the series of transitions between secondary structures. This series of transitions or microscopic structure can be described by a system of ordinary differential equations that can be approximating using balanced truncation to determine the macroscopic structure. By doing so, we have been able to identify signature topological features of microscopic structure and mathematically characterize the corresponding classes of macroscopic structure. Thus we are now able to take large, complex systems, reduce them, and understand their behavior. In the future, we hope to be able to identify small microscopic changes that lead to large macroscopic changes and possibly phasetransition like conditions between secondary states. Ultimately, this may lead to the development of a secondary-structure kinetics theory describing how one or more strands of DNA pair with one another to form different secondary structures and its potential future experimental verification. Item Type: Thesis (Senior thesis (Major)) Subject Keywords: (Physics) Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Physics Thesis Availability: Public (worldwide access) Research Advisor(s): Mabuchi, Hideo Thesis Committee: None, None Defense Date: 8 May 2006 Record Number: CaltechTHESIS:02132018-142852598 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:02132018-142852598 DOI: 10.7907/4bxt-wa89 ORCID: Author ORCID Meruelo, Alejandro Daniel 0000-0001-6087-1818 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 10718 Collection: CaltechTHESIS Deposited By: Joy Painter Deposited On: 13 Feb 2018 23:49 Last Modified: 14 Nov 2025 19:47 Thesis Files Preview PDF - Final Version See Usage Policy. 19MB Repository Staff Only: item control page

From Sequence to Function through Secondary Structure Kinetics of RNA and DNA Thesis by Alejandro Meruelo In Partial Fulfillment of the Requirements for the Degree of Bachelor of Science California Institute of Technology Pasadena, California 2006 (Defended May 8, 2006) '-l 1 } ' • 11 © 2006 Alejandro Meruelo All Rights Reserved '( . iii Acknowledgements I would like to thank my mentor, Professor Mabuchi for his guidance in selecting reading relevant to my research, for fruitful discussion, and for giving me the opportunity to pursue this project. I would also like to thank Ben Rahn for the seed materials for the project, and Kevin McHale and Nathan Rodas for very helpful discussion; Professor Mabuchi and Kevin McHale were also very helpful in reviewing this document. This project was made possible by the support of the Mellon Mays Undergraduate Fellowship (MMUF). IV Abstract A number of efforts to determine function from sequence of RNA and DNA have been made with varying success. Here we study the determination of function from sequence of DNA and RNA through their secondary structure kinetics, specifically the series of transitions between secondary structures. This series of transitions or microscopic structure can be described by a system of ordinary differential equations that can be approximating using balanced truncation to determine the macroscopic structure. By doing so, we have been able to identify signature topological features of microscopic structure and mathematically characterize the corresponding classes of macroscopic structure. Thus we are now able to take large, complex systems, reduce them, and understand their behavior. In the future, we hope to be able to identify small microscopic changes that lead to large macroscopic changes and possibly phasetransition like conditions between secondary states. Ultimately, this may lead to the development of a secondary-structure kinetics theory describing how one or more strands of DNA pair with one another to form different secondary structures and its potential future experimental verification. t ' v Contents Acknowledgements iii Abstract iv 1 Introduction 1 1.1 Problem Statement . . . . . 1 1.2 An Input-Output Approach 4 2 The Model Reduction Algorithm 2.1 Principles 2.2 Speed 5 8 . . 12 2.3 Implementation 14 2.4 The Statistics of Model Reduction . 18 3 The Mathematical Characterization of Macroscopic Structure 4 8 21 3.1 3-state Systems 21 3.2 4-state Systems 31 3.3 The Trapping-state System . 37 3.4 Distance between Two Digraphs . 40 The Signature Topological Features of Microscopic Structure 42 4.1 The Loop . . . . 42 4. 2 The Linear Chain 43 Discussion 45 ". VI 5.1 Theoretical Implications 45 5.2 Future Work . . . . . . . 46 Bibliography 49 Appendices 50 A Code for Model Reduction 51 f I Vll List of Figures 1.1 Configmation space of the DNA strand AATT, excluding pseudoknots. 2 1.2 Digraph for the DNA strand AATT. 2 1.3 The big picture. 4 1.4 A simple 3-state system. 6 1.5 The model reduction scheme. 7 2.1 The observability and controllability ellipsoids. 10 2.2 The balancing transformation: alignment of the observability and controllability ellipsoids. . . . . . . . . . . . . . . 11 2.3 Plot of t(s) v. n for the balancmr algorithm . . 13 2.4 Block diagram of steady-state removal and restoration procedure. 15 2.5 Model reduction distribution under limited error bound. 19 3.1 First 3-state system. 25 3.2 Second 3-state system. 27 3.3 Third 3-state system . . 27 3.4 Fourth 3-state system. 28 3.5 Fifth 3-state system. 28 3.6 Plot of r12 = r12(r21, r32). 29 3.7 First 4-state system. 32 3.8 Second 4-state system. 33 3.9 Third 4-state system . . 33 3.10 Fourth 4-state system. 34 3.11 Fifth 4-state system. 34 ! ' Vlll 3.12 The trapping-state system. 37 3.13 Parameter-space of the trapping-state system. 38 3.14 Bounded area of the trapping-state system. 39 4.1 Reduction of a loop. . . . . 43 4.2 Reduction of a linear chain. 43 5.1 Example of a large system. 45 •• 1 Chapter 1 Introduction 1.1 Problem Statement This work during the past summer and academic year has concerned understanding the secondary structure kinetics of RNA and DNA. Understanding protein folding for a long time has posed a problem of a similar nature to understanding the secondary structure kinetics of RNA and DNA. Theoretical approaches to understanding protein folding have centered around the development of models of free energy landscapes containing the different intermediates through which a protein travels to achieve its lowest-energy conformation. These models have proved crucial to protein structure prediction and design. In the same way that theoretical approaches have been taken to understanding protein folding, we discuss a novel theoretical approach for understanding the secondary structure kinetics of RNA and DNA and its potential experimental implications. Understanding DNA and RNA hybridization kinetics at the level of secondary structure has promising potential application to molecular dynamics and nano-biotechnology. For example, in molecular biology, it may help us better detect pathogens, mutations, and play an important role in drug screening [l ]. It may also be crucial to understanding the activity of DNA and RNA catalysts [2]. Finally, this understanding will have important consequences for nanotechnology and DNA-based computing, e.g. for designing devices like DNA tweezers [3]. A number of people have set the stage for important work in this area. Some have .... 2 compiled databases of thermodynamic parameters for many motifs for the purpose of making accurate DNA structure prediction [4] and have implemented programs capable of prediction of secondary structure for pairs of RNA or DNA strands using these data like Pair Fold (http://www.rnasoft.ca/cgi-bin/RNAsoft/PairFold/pairfold.pl). Our work differs from these in that we hope to be able to understand the theory underlying secondary structure kinetics and thus to be able to draw larger conclusions about the way that it occurs. A- A- T- T T- T ! I A- A T I T I! T- T (T:;lf\ A- A f- 0J A Figure 1.1: Configuration space of the DNA strand AATT, excluding pseudoknots. A- A- T- T Figure 1.2: Digraph for the DNA strand AATT. We begin by considering one or more strands of DNA and the fashion in which they hybridize with one another. Consider the DNA strand AATT, for example. We will be considering all possible base-pairings, excluding pseudoknots (Figure 1.1). An example of a pseudoknot structure has been excluded in the figure shown. We assume that hybridization can occur only through certain elementary kinetic •• 3 steps. These are (i) the formation of a new base-pair, (ii) the breaking of an existing base-pair, and (iii) a shift from one base pair to another. Then we may depict the configuration space, containing all possible configurations of a set of DNA strands, as a graph with directed edges (Figure 1.2). Each node of the graph then corresponds to a possible secondary structure and each single-step transition is determined by thermodynamic considerations. Secondary structure dynamics can consequently be viewed as a random walk on this graph. The central question of my research and the question I will address is, "How do we take this graph and simplify it?" The goals of our research are to apply a model reduction scheme to this graph in order to arrive at a more intuitive understanding of large systems (with associated microscopic structure) so that we can predict the way hybridization occurs. Specifically, we aim to be able to (i) identify topological features of microscopic structure that correspond to particular classes of macroscopic structure (of reduced systems) and (ii) to be able to mathematically characterize these. The larger goal will be to identify small changes in microscopic structure that lead to substantial changes in macroscopic structure. By gaining a greater grasp of the manner in which hybridization occurs, we then hope to be able to learn more about things governed by secondary structure kinetics. We can begin by answering questions concerning, e.g. the time dynamics of ribozymes, RNA-based enzymes. We aim to be able to understand the mechanisms that stabilize specific structural folds. Furthermore, from studying the manner in which hybridization progresses, we hopefully will be able to create new RNA and DNA-based structures with specific functional behavior for the purposes of bioengineering, e.g. for gene therapy. In addition, we hope that the theory we develop will help us understand disease processes involving pathogenic ribozymes in plants and animals such as Hepatitis delta virus. Concretely, we will be determining how long and with what likelihood it takes to transition from one secondary structure to another, with the ultimate goal of being able to predict the manner in which hybridization progresses starting with different configurations of DNA or RNA strands. We lay the foundation to accomplish the . ' 4 goals we have outlined here in the next few chapters. 1.2 An Input-Output Approach macroscopic picture hltp:l/chemistry.ucsc.edu/-wgscott/ physical system microscopic picture system dynamics Figure 1.3: The big picture. We propose approaching the problem of understanding the secondary structure kinetics of RNA and DNA in the following fashion. Imagine a physical system, such as that depicted in the initial sequence of Figure 1.3, whose time dynamics we would like to better understand. We will be interested in studying the corresponding microscopic picture underlying the transitions required to go from one secondary structure to another for this physical system. We can depict them using a directed graph where each node represents one possible secondary structure and the arrows connecting them the possible transitions between these structures (Figure 1.3). We may then approximate this microscopic picture by a simpler macroscopic picture, as shown in the figure to the right. Through such an approach, we hope to be able to gain insight into how hybridization occurs. More concretely and differently than others have done in the past, we will be treating our physical system as an input-output system. We are interested in studying the time dynamics for going from one initial state to some final state. So, the input will be the initial secondary state of the system and the output will be the final secondary state of the system. .. 5 We begin by considering a system with time varying input u(t) that drives a system with state x( t) and output signal y( t). We examine the input-output map and the system dynamics, which can be written as ±(t) = J(x(t), u(t)) (1.1) y(t) = h(x(t), u(t)). We model the forward and reverse transitions between individual secondary structures each by a single rate, defining the transition-rate matrix A in which the coefficient rii defines the transition rate from state i to j: ;tx = Ax(t), -I:rj1 r12 r1n -I:rj2 #2 r2n rn2 -2.:rjn j=Jn jfl r21 A= rn1 (1.2) In general, the matrix A will not be symmetric. Note that all rates are by definition nonnegative. Because our forward and reverse transitions are defined by a single rate, we may then write our system in matrix form as the linear system: x = Ax+Bu, y=Cx, (1.3) where A,B, and Care time-independent real matrices of sizes n x n, n x m , and p x n, respectively. B represents the initial secondary state (input of the system) encoded in the form of an initial probability distribution vector ... 6 1 0 B= 0 corresponding to the first state and C represents final secondary state (output of the system) encoded in the form of a final probability distribution vector C=[o O... 1] corresponding to the last state. Figure 1.4: A simple 3-state system. To make this concrete, consider the example of the simple 3-state system depicted in Figure 1.4. We can write down its corresponding matrices as A= 0 ri2 0 -(r12 + r32) 0 r32 0 1 ' B= 0 0 Next, to solve the approximate time dynamics of traveling from one initial state to some final state for a large system representing our model of single-step DNA hybridization of order n, we will use two input-output model reduction algorithms from control theory. We illustrate these in Figure 1.5. We will be applying these two algorithms to arrive at a lower order approximate model of the behavior of our system. The first of these will be the Kalman decomposition, which eliminates any degrees of freedom that don't connect the initial to the • I ·- I 7 Kalman decomposition (no loss of accuracy) balanced truncation (bounded error~ ~1 1c - c" - 11<_, ~2IoJ ::::: ~ Figure 1.5: The model reduction scheme. final state to arrive at an n' order model. The second will be balanced truncation, which essentially rank orders the degrees of freedom of the problem and truncates the degrees of freedom that have the least effect on input-output characteristics to arrive at a simpler model of order n". This is possible because there often will be a wide gap between the least and most important degrees of freedom, allowing for the identification and removal of those less contributing degrees of freedom. The latter algorithm has a complexity of O(n 2 ) in space and O(n 3 ) in time. The mathematics underlying these two algorithms is discussed in some depth in chapter 3. It turns out that the model reduction scheme we are performing has a bounded error as we will see in the next chapter, giving us great control over our ability to capture the dynamics of our original model with a well-defined degree of accuracy. Our model reduction technique as a whole is particularly important because large systems are unwieldy to solve and understand, while small ones are often not. By being able to reduce large systems, we will be able to gain intuition into how they behave. We hope to be able to understand how particular microscopic topological features reduce to certain classes of macroscopic structure and then to be able to mathematically characterize how such systems evolve in time. We hope that our approach in the long run will be adapted to understand much larger scale systems and more complicated biological systems. 8 Chapter 2 The Model Reduction Algorithm In this chapter, we will review in some depth the model reduction scheme that we will be applying to understand and make predictions about the manner in which hybridization occurs. Specifically, we will address the questions, "How does the model reduction algorithm work?", "How fast does it work?", "How did we implement it?", and "How well does it reduce the systems we will be encountering?" By addressing these questions, we hope to be able to provide some insight into the underlying methods of our research and to motivate the chapter to follow. 2.1 Principles We begin by going over the mathematical principles underlying our model reduction algorithm. As briefly touched upon in Chapter 2, we will be applying two algorithms in succession to our system, the Kalman decomposition and balanced truncation. Our methodology is based on the work of Rahn [5], a former graduate student in the lab. We first discuss how the Kalman decomposition algorithm works. Any dynamical system can be described by a variety of coordinate systems (e.g. a simple pendulum can be represented in terms of the angle of the string with the vertical or its vertical displacement). The Kalman decomposition reduces a given system to its lowest order, exact model in order to simplify system analysis. We begin by examining the input-output map '11 : u r-+ y neglecting any internal behavior. Then we introduce the concepts of observability and controllability. If a . ' 9 state is unobservable, it does not affect the output and if a state is uncontrollable, it is unaffected by the input. A minimal realization is a model in which all states are both controllable and observable. The goal of the Kalman decomposition is to obtain a minimal realization. To do so, we will need to equate our system to the controllable and observable part of our dynamical system. The Kalman Canonical Decomposition Theorem [6] states that given a dynamical system described by equation (1.3), a similarity transformation T exists such that the equivalent linear system [TAT~' cr- 1 T: j - Aco 0 A21 A:a .A23 .A24 Eco 0 0 Aeo 0 0 0 0 .i143 Aco 0 Ceo 0 Ceo 0 .i113 0 0 Eco can be computed, where c stands for controllable and o for observable (and their complements c,o for uncontrollable and unobservable, respectively). More importantly, however, we are interested in the transfer matrix, i.e. just the input-output behavior with no regard for internal behavior. Thus we have that the transfer matrix of the system is equal to its controllable and observable part, i.e. A B [ C 0 J _ [ ~co Eco Ceo O j Hence we have our minimal realization and the lowest order, exact model for our system. Next, we will look at how the balanced truncation algorithm works. The exact balanced truncation algorithm is well-suited to problems of small to medium size and has complexity O(n2 ) in space and O(n 3 ) in time [8]. Intuitively, balanced truncation ._' 10 involves identifying those degrees of freedom that are only weakly changed by the input and most weakly affect the output, and removing these. It is often handy when building simulations and theoretical models for the evolution of macroscopic quantities in large, complex systems, as we have here. We begin performing balanced truncation by first introducing the observability gramian of (C,A), called Yo, and the controllability gramian of (A,B), called Xe, which are matrices that provide useful quantitative measures of a state's observability and controllability: 1 1 00 Yo = eA*r C *CeAr dT 00 Xe= eAr BB*eA*r dT These terms constitute the solutions of the two Lyapunov equations [6) A*Ya + YaA + C*C = 0, AXc + XcA* +BB*= o. Figure 2.1: The observability and controllability ellipsoids. Alternatively to using these definitions and more intuitively, we may geometrically interpret the observability and controllability gramians by the introduction of the observability and controllability ellipsoids (Figure 2 .1), i.e. 1 Ea= {Yo2 xa: Xa E Cnandlxal = 1} 1 Ee= {XJxc: Xe E Cnandjxcl = l}. . ' 11 Figure 2.2: The balancing transformation: alignment of the observability and controllability ellipsoids. Then the largest prinicipal axes lie along the most observable or controllable directions. In order to arrive at a good approximate model, it would make most sense to eliminate the least observable and controllable states since these are the least important. By performing a balancing transformation (requiring A stable, i.e. Re(..\(A)) < 0), which results in the simultaneous diagonalization of the observability and controllability gramians, we align the observability and controllability ellipsoids (Figure 2.2), thus making sensible truncation possible. The diagonalization of both gramians under congruence transformations is made possible by the fact that the gramians are Hermitian and positive [7]. The resulting diagonal matrix is of the form 0'1 Thus we have a diagonal matrix with decreasing Hankel Singular Values (HSVs) along the diagonal. In order to arrive at an approximate model of order n", we must truncate the system only where the HSVs have a clear separation (i.e. are distinct), beginning with the smallest of such states, which are the least controllable .. 12 and observable. Partitioning our matrix as mn" x n" A n E i& our reduced model becomes This completes the process of performing balanced truncation on our system and the approximate n" order model produced preserves the stability of our system with bounded error k 0-1 ::=; llG - G~11 Iloo ::=; 2 L O"i, i=l where the o-i correspond to the distinct truncated HSVs [5]. Thus, by understanding the underlying mathematics of these two algorithms, we can better understand why they are so effective for our purposes, as well as have a greater appreciation for how the approximation scheme we have employed works. 2.2 Speed The next most pressing concern is, "How fast can we reduce these types of systems?" This will give us an idea of how large a system we can handle and how practical our approach is for tackling these types of systems. 13 Data for the balancmr algorithm (N=lOO) 0 t (s) n 0.540 20 1.230 50 16.230 174 231.910 410 100 200 3.00 I 400 X=n Figure 2.3: Plot of t(s) v. n for the balancmr algorithm. We benchmarked the longest step of our model reduction scheme balancmr, a cornrnand in MATLAB (Mathworks Inc., Natick, Massachusetts, USA) that performs balanced truncation, to see how large a system it could reduce to a 3-state system in times of order 0.5s, ls, 10s, and lOOs. We chose to use the rss (n, p, m) function provided in MATLAB to generate a random, stable, nth order model with m inputs and p outputs. Specifically, we ran rss (n, 1, 1) to generate a random stable model of order n with one input and one output. Running N=lOO trials and then using the profile command on a 1.5 GHz Power PC G4 with 1 GB of RAM, running Mac .. 14 OS X 10.4.3 and MATLAB 7.0.4.352 (R14) SP2 we obtained the results given in the table above and in Figure 2.3. From the figure, we conclude that with increasing time, the pace at which higher order models can be reduced using balanced truncation increases. One question that we raise for experimentalists to answer is, how do our randomly generated systems compare to typical systems involved in studying the secondary structure kinetics of DNA or RNA and how can we simulate more realistic systems? The study of the speed of model reduction as a whole is an interesting one and one which underlines the important practical constraints of applying this strategy to answer research questions in general. Since we are currently in the process of developing this technique, we will typically be dealing with systems of order 50 at maximum and thus conclude that our model reduction scheme is sufficiently fast for our purposes. In the future, it is valuable to know that the algorithms we are using can be adapted. For example, the balanced truncation algorithm, while in its exact form is well suited for small to medium problems like our current one, can be adapted to larger problems by approximate application [8]. This will be especially important since the size of DNA configuration space grows exponentially. 2.3 Implementation I now discuss the implementation of the model reduction algorithm examined in the first section of this chapter. The corresponding MATLAB code may be found in the Appendix. Specifically, I will comment on how the two algorithms for reduction were implemented and what hurdles had to be overcome in order to accomplish the goals of the model reduction scheme as a whole. The implementation of the two algorithms, the Kalman decomposition and balanced truncation, is a straightforward process in MATLAB, which is equipped with toolboxes containing these functions. These come in the form of a minimal option for constructing a state-space model using the ss (stable_matrix, 'minimal') command and the balancmr(stable_matrix,3) command to reduce a given system to .. "' 15 Figure 2.4: Block diagram of steady-state removal and restoration procedure. one of order 3, respectively. However simple these may seem, the balanced truncation algorithm requires that the transition-rate matrix A be stable, i.e. Re(.A(A)) < 0. This means that there is a need for the removal of the steady state that arises as a consequence of the physical model of our system and otherwise invalidates the stability condition. Also, we will need to to restore the system once it has been reduced using the balanced truncation algorithm. This procedure as a whole has been illustrated in Figure 2.4. We begin by outlining the mathematical problem we will solve. Consider the linear state-space input-output model of earlier: x = Ax+Bu, y=Cx, (2.1) By the nature of our problem, A is a matrix containing one zero eigenvalue, corresponding to the steady-state and (n-1) negative eigenvalues, corresponding to the transient states. Performing balanced truncation on our system naturally requires that A be stable, i.e. Re(.A(A)) < 0. Then we must develop a method to remove this steady-state value and later restore it after balanced truncation has been performed. Our first step will be to diagonalize - or in the case of degeneracy, find the Jordan form J of our matrix A. For simplicity, assume that all our eigenvalues are distinct so ~ • • t 16 that we will be diagonalizing A using a transformation T: D = r- AT 1 ) where D is a diagonal matrix made up of the eigenvalues of A, and T is a matrix whose columns correspond to the eigenvectors of A and are ru:ranged in an order corresponding to their eigenvalues appearing in the matrix D. Suppose that the zero eigenvalue appears in the qth row and column of D. Then we can define a new matrix V as the rectangular n x (n - 1) matrix resulting from deleting the qth column of T and v- 1 as the rectangular (n - 1) x n matrix resulting from deleting the qth row of r- 1 . Then we have a transformation v- 1 that maps from the original state space onto the stable subspace of A. To apply this change of basis we substitute i (and x =Vi) into equation (2.1) to obtain i =: =v- ± =v- (Ax +Bu)= v- Avx + v- Bu =Ax+ Bu 1 1 1 1 y = Cx = CVx =Ci, which then yields that A= v- 1 AV,B = v- 1 B,C =CV Thus we have succeeded in obtaining a transformation to remove the steady-state dimension of the state-space. Once we've performed balanced truncation on our system to obtain Xr, we are then ready to restore our system to its original form by reincorporating the steadystate value. To do so, we'll use an expression similar to (2.1) above for our reduced system: ±o = IIoBu, Y = CrXr + Cxo, u = o(t) 4 • • .. 17 where I1 0 is the orthogonal projector onto the steady-state (zero eigenvector) of A. Then fort> 0, since only these times make physical sense, xo = IloB. To compute I1 0 , it may be easier to think back to quantum mechanics where we could write B = L li)(i'jB. i Then we may write B = (Ilt + Ilo)B = IltB + I10 B, where t stands for transient and 0 refers to the steady-state. This simply means that Ilt = L li)(i'I and Ilo = !O)(O'j. i>O Then we can now write our final expression for the output Thus by introducing the factor CI10 B, we are able to restore the steady-state contribution to our reduced system. We now discuss the code developed to formally accomplish this procedure. To achieve the goals discussed above, the zero eigenvalue of the transition-rate matrix A, corresponding to the steady-state of the system, must be removed. The code included in the appendix begins by looking for this zero eigenvalue, constructing an orthogonal projector onto the steady-state (zero eigenvector) that we will need later on, making a map from the original state space to the stable subspace (by throwing out the components corresponding to the steady-state), and then performing the transformation from the original state-space to a strictly stable one. It then applies the Kalman decomposition to the resulting matrix using the ss (stable_matrix, 'minimal') command. .. ' Ill .. 18 This is followed by performing balanced truncation on the resulting minimal order, exact model, reducing it to a system of order 3. It is at this point that the steadystate behavior is restored using the orthogonal projector constructed earlier and the system is once again a good approximation of its original functional behavior. The steady-state removal and restoration hurdle was a major problem that had to be overcome and there still remains the need to restore the reduced system to continuous-time stochastic (CTS) form. The general idea behind this problem is given a reduced matrix A, "Can we restore it to CTS form (for which each column sums to zero, and typically diagonal elements are negative and off-diagonal elements are positive or zero) while simultaneously preserving its current eigenvalues?" This is a difficult problem not solved by former graduate student Ben Rahn and is currently being worked on by another graduate student, Kevin McHale. The most that can be said right now about this problem is that for an n x n CTS matrix, all of the eigenvalues must satisfy the relation llm(>.)J < t(7r) IRe(>.)I - co n. The proof of this relation is due to Bryan Shader. One idea I had to solve this problem would involve Lagrange multipliers through which the problem's constraints could be incorporated. However, as a whole, we have succeeded in producing a system whose functional behavior approximates the original, though it has lost the structural features of its digraph by coming out of CTS form. Thus we are able to compute effective transition-rates, but not yet study the digraph structure of these systems with confidence. 2.4 The Statistics of Model Reduction It is appropriate to close this chapter by motivating the following chapter. We are interested in looking at the way large systems reduce to smaller systems. Specifically, we will be interested in the typical reduction of order under a fixed error bound. Recall that the error bound for systems using our model reduction scheme is given . ' 19 by k 0"1 S jjG - Gn lloo S 2 11 L O"i, i=l where the ai are the distinct truncated Hankel Singular Values (HSVs) of the system [5]. We will be interested in what typical order models we obtain as a result of reducing our system when the error is approximately 20% of the system as a whole. 8 7 6 5 4 3 2 1 0 I- Q) "O I- Q) I- Q) I- "O 0 I- "O 0 0 I- "- Q) "O "- 0 • Initial Order 20 • Initial Order 30 I- Q) "O I- 0 I'-I- Q) "O I- 0 I- Q) "O "- 0 00 I- Q) "O "- 0 I- Q) "O "- 0 ('I) I- Q) "O I- 0 I- Q) "O I- 0 ..-1 I- Q) I- Q) "O "O 0 0 I- I- Figure 2.5: Model reduction distribution under limited error bound. In general, we will at worst be dealing with a system of order 50 and more likely will be dealing with a system of order 20 or 30. I analyzed the distributions of model reduction of systems of order 20 and 30, with a sample size of N ;:::::: 20 for each. To accomplish the goals of my study I used the MATLAB command rss (30, 1, 1), e.g. to randomly generate a 30-order model, the command minreal to reduce it to its lowest order exact model (Kalman decomposition), and the command balancmr to truncate it according to its HSV plot (balanced truncation). The criterion for truncation was that the sum of the truncated HSVs was closest to 103 of the total sum of HSVs. The resulting distribution plots for initial models of order 20 and 30 are given in Figure 2.5. The associated statistics for these plots have been tabulated in the table below. 20 Model reduction distribution statistics Initial Order (N) 20 (21) 30 (22) Average 5.238 3.955 Median 5 4 Standard Deviation 2.256 1.362 From these results, we can infer that typically under the 203 error bound we will find that reduced models will be of order 4 or 5. Whether our systems will remain of this order when restored to CTS form is another question that must be answered in the future. Presuming this is so for simplicity of discussion, since many order 5 and some order 4 systems cannot be solved analytically, we will begin by analyzing all possible 3-state systems and then learn as much as we can from those 4-state systems that can be solved analytically in order to mathematically characterize the behavior of typical systems we will be obtaining using our model reduction scheme. 21 Chapter 3 The Mathematical Characterization of Macroscopic Structure 3.1 3-state Systems We begin by analyzing 3-state systems with the intent of gaining intuition into larger systems that reduce to such systems using our model reduction scheme. Specifically, we will be computing the probability P1i(t) of the transition from the initial state i to final state f, where To solve for P1i(t), we will have to solve the problem set up in Chapter 2 consisting of the system of equations (1.3), where P1i(t) will simply be the output of the system y( t). We restate these equations here along with some additional conditions for our problem: x=Ax+Bu, y=Cx, u = o(t) , x(O) = Xo , x E ~n, u E ~m, y E ~P We would like to solve this problem for t ~ 0 since no other times have physical meaning. We begin by solving the homogeneous equation 22 i: =Ax whose solution is well known to be Next we look at the inhomogeneous equation i: = Bo(t) whose solution is x(t)i = BB(t), where 8( t) is the Heaviside function defined by 1, t ~ 0, B(t) - { 0,t < 0. Since we are interested in time t ~ 0 we will simply take our solution to be x(t)i = B, t ~ 0. Adding our two solutions to obtain the general solution we have x(t) = eAtx0 + B + C. Since x(O) = x 0 then x(O) =x0 +B+C;C= -B so that our general solution becomes However, we know that at time t=O, the initial state is specified by the initial probability distribution vector B so that x(t) = eAt B, ; y(t) = CeAt B = PJi(t), t ~ 0. 23 Thus we have found the solution for P1i(t) that we will use for the remainder of our work. Exponentiation of the matrix A is relatively straightforward to accomplish by solving for the corresponding eigenvectors, from which one can construct a basis by which the original matrix can be diagonalized by its eigenvalues and then easily exponentiated and transformed back using the original basis. In the case of degeneracy, this task can be accomplished using the Jordan canonical form, which also gives insight into the form of the solution. By mathematically characterizing the functional behavior of these simple systems, we will learn about their characteristic behaviors and thus will be primed to understand systems that reduce to these simpler systems when we encounter them later on. Before we introduce the main results of this chapter, it is useful to note some of the more general properties of n-state systems that we have learned by studying 3- and 4state systems. We have found that the functional form for the time dependence of the transition from an initial to a final state has a number of characteristic properties. The functional form is dependent solely on the number of distinct or degenerate eigenvalues making up the system and their degree of degeneracy. The number of secondary states of a system determine the number of eigenvalues for a given system. It will be useful to term the eigenvalues of the system its eigenrates. Then these eigenrates are given by expressions involving the transition rates making up the system itself. Note that the eigenrates have been generally found to be real for all systems studied. When all of the eigenrates for a system are distinct (and real), the functional form for its time dependence consists of a linear combination of one exponential term for each eigenrate ->.i of the form e->.it. When an eigenrate is degenerate (and real), any number of similar exponential terms multiplied by a polynomial of increasing order may be added to the first purely exponential term up to a maximum order of one less than the order of degeneracy. This fact is a corollary of the Jordan Canonical Form Theorem, but can also be understood by inspecting the general form of an exponentiated matrix of Jordan form J multiplied by time, eJt [9]. To consider an 24 example, if we have a 5-state system with a degeneracy of 3, the functional form for the transition rate P51 ( t) from the initial to the final state will be where the eigenvalue 0 always occurs and is ignored in the counting of eigenvalues for the system. Each of the coefficients {a, b, c, d} is often a complex expression involving the transition rates of the system and is determined according to the structure of the digraph (e.g. Figure 3.1). Some of these may vanish, allowing for more flexibility in the number of functional forms resulting from the same digraph structure. We now introduce some of the main results of the chapter. We begin by characterizing the functional behavior of all possible 3-state systems. In constructing these systems, we assume that all systems have a single absorbing final state. This means that once in the final state for a system, there is zero likelihood of exiting this state to return to any other state in the system. Here we enumerate all observed classes of functional behavior for a 3-state system resulting from analytically solving for their time dynamics. These classes of behavior were discovered by using rate-limiting conditions (e.g. .:\ 1 » .:\2 ) and degeneracy conditions (e.g . .:\ 1 = .:\ 2 ) to produce different functional classes across many different types of systems studied. Note that these have been written using the minimum number of parameters based on the fact that P 31 (t = 0) = 0: Possible classes of functional behavior (3-state systems) 1 - (1 + b)e-:>.it + be- >-2 t 1 - e->.it + be->- 1 tt Note the resonance (involving a polynomial) class of behavior that results from degenerate eigenvalues as distinguished from the other class of behavior that results from distinct eigenvalues. Next, we examine the 5 possible 3-state systems under the assumptions above. Before going on to state what we observed about these systems, we give the eigenvalues .:\ 1 , .:\ 2 and the coefficients for these systems, including two examples of how these are 25 determined. For every system, we will provide two solutions (or state that one doesn't exist), one for the case of distinct eigenvalues with general form 1-(1 +b)e->- 1 t+be->- 2 t and one for the case of degenerate eigenvalues with general form 1 - e->-it + be->- 1 tt. Note that some systems have no degenerate solution for positive rates. Figure 3.1: First 3-state system. Beginning with the system depicted in Figure 3.1 , we write down its corresponding matrices given by 0 A= 1 0 ' B= ,c=[o o 1]. 0 0 0 The eigenvalues of our system are visibly given by >. 1 = r 21 and >. 2 = r 32 . We start by considering the distinct case (i.e. >. 1 =I >. 2 ). Recall that P31 (t) = CeAt B. For the distinct case, we will exploit the ability of the matrix A to be diagonalized by its eigenvalues >.i and then easily exponentiated. Specifically, eAt = Sdiag[e->-it]s- 1 , where S is the similarity transformation yielding the diagonal form of the matrix A. Constructing this similarity transformation using the eigenvectors of the matrix A we find: 0 -l+!ll. r 32 S= 0 _!ll. 1 1 r 32 0 -1 1 , diag[e- >.it] = 1 0 0 0 e-r21t 0 0 0 ' 5-1 e-r32t Then we may find P 31 ( t) by multiplying the above matrices: = 1 1 1 ---2:.ll.._ 0 0 -1 0 r21-r32 _---2:ll_ r21-r32 26 Comparing this to the general distinct form expected 1 - (1 + b)e->.it + be->- 2 t we state the parameters of our solution for this system to be: Ai = r2i ' A2 = r32 ' b = --2:ll_ r32-r21 Next, we consider the degenerate case (i.e. Ai = A2 ) for the same system with corresponding matrices given above and eigenvalues Ai = r 2 i and A2 = r 32 . We begin by setting r 32 = r 2 i in the matrix A. As before, P3 i(t) = CeAtB. Since we are dealing with degeneracy, we will employ the Jordan form J. Then we may write eAt = S eJt s-i, where Sis a similarity transformation yielding the Jordan form. Using an algorithm for finding a similarity transformation S constructed from a basis of eigenvectors such that a degenerate matrix A with real eigenvalue A= r 2 i assumes its Jordan form we find [9]: 0 S= 0 0 -1 1 1 1 r21 i r21 0 ' ]= 0 0 0 0 -T2i 1 0 0 -r2i ' s-i = 1 1 1 -1 -1 0 -T21 0 0 Using the well-known structure of a Jordan form when exponentiated we obtain: 1 0 0 0 0 Thus we are now prepared to compute the value of P 31 (t) by multiplying the above matrices and find: 27 Comparing this to the general degenerate form expected 1 - e-.\it + be-.\ 1 tt we state the parameters of our solution for this system to be: Figure 3.2: Second 3-state system. Using this method again, for the next system in Figure 3.2 we found that its distinct solution is given by the parameters: >.1,2 = !(r12 + T21 + T32 ± .J(r12 + T21 + T32) 2 - 4r21T32), b= r r r21 r32-(r12+r21 +r32)(r12+r21 +r32- (r12+r21 +r32) 2 -4r21 r32)+ ~(r12+r21 +r32-y' (r12+r21 +r32) 2 -4r21 r32) 2 A degenerate solution for this system is not possible due to algebraic constraints (an example of such a constraint is given later in this chapter). Figure 3.3: Third 3-state system. For the system in Figure 3.3 we found that its distinct solution is given by the parameters: >.1,2 = !(r12 + T21 + T31 ± .J(r12 + T21 + T31) 2 - 4r12T31), r31(r12-r21-r31+y'(r12+r21+r31) 2 -4r12r31) b= 2(r12r31 +(r12+r21 +r31)(r31-r12-r21 +y'(r12+r21 +r31) 2 -4r12r3i)+~(r12+r21 +r31-y'(r12+r21 +r31) 2 -4r12r31) 2) A degenerate solution for this system is not possible due to algebraic constraints. For the fourth system shown in Figure 3.4 we found that its distinct solution is given by the parameters: Its degenerate solution is given by the parameters: 28 Figure 3.4: Fourth 3-state system. Figure 3.5: Fifth 3-state system. Finally, for our last 3-state system we found that its distinct solution is given by the parameters: >-1,2 = Hr12+r21 +r31 +r32±J(r12 + r21 + r31 + r32) 2 - 4(r12r31 + r21r32 + r31T32) ), b _ r12r31 -r21r31-r~ 1 +2r21r32+r31r32+r31 a (r12+r21 +r31 +r32) 2 -4(r12r31 +r21r32+r31r32) 2(o+.6+'Y) ' = T12T31 + r21 T32 + T31 T32 /3 = (r12+r21 +r31 +r32)(-(r12+r21 +r31 +r32)+J(r12 + r21 + r31 + r32) 2 - 4(r12r31 + r21r32 + r311 'Y = ~(-(r12+r21 +r31 +r32)+J(r12 + r21 + T31 + r32) 2 - 4(r12T31 + r21r32 + T31T32))2. A degenerate solution for this system is not possible due to algebraic constraints. Next we look at the number of free parameters and classes of functional behavior observed for these systems. The number of classes of functional behavior was determined by studying rate-limiting conditions (e.g. >.. 1 » >.. 2 ) and degeneracy conditions (e.g. >.. 1 = >.. 2 ), and then numerically simulating that the associated classes could indeed be found for each of the digraphs. 29 To give an example of a very simple rate-limiting condition, consider the fourth system (Figure 3.4). Suppose that we would like to simulate the rate-limiting condition >-2 » ..\ 1 . Looking at its eigenvalues, this simply means that r 32 » (r21 + r 31 ). Numerically setting r 32 = 107 and r 21 = r 31 = 1 we obtain So we've been able to induce the class 1 - (1 + b)e->- 1 t + be->- 2 t, with b=O. Ratelimiting conditions can become significantly involved as one can see from looking at the eigenvalues given above for some systems. More specifically for more complicated systems, one must solve for an algebraic condition involving one rate compared to the others and then numerically check that simulating such a condition results in the expected outcome. Sometimes when one rate-limiting condition fails, another may still work to induce a particular class. Inducing degeneracy is performed similarly by solving for a particular rate in terms of the other rates and can prove tricky as well. Sometimes a condition can prevent degeneracy from ever occurring. To give an example of a simple failure to produce degeneracy, consider the second system (Figure 3.2). We would like to induce the degeneracy condition ,\ 1 = ..\ 2 . Looking at the eigenvalues, this means that (r 12 + r 21 + r 32 ) 2 = 4r21 r 32 or r 12 = -(r21 +r32 )±2Jr21 r 32 . We hypothesize that the rate r 12 is never positive (assuming the rates r 21 and r 32 are positive as required by our model) by observing that this is the case in Figure 3.6 in which the two possible functions for the rate r 12 are plotted as a function of r 21 and r 32 • To show that this is true, let us maximize r 12 to see 30 whether we can obtain a positive value. Critical points occur at r 21 = r 32 and r 21 = 0 or r 32 = 0. But we require that r 21 > 0, r 32 > 0, eliminating the latter two points. Then at r21 = r32, r12 = -(r21 + r21) ± 2r21 = 0 or < 0. But r 12 > 0. Thus ..\ 1 i= ..\ 2 and the corresponding degenerate class 1 - e->-it + be->- 1 tt is precluded by means of this algebraic constraint. The work outlined in this chapter as a whole involved solving a total of 10 systems analytically for their associated eigenrates made up of all the different transition rates for each system and then numerically simulating over 80 rate-limiting conditions to verify whether they led to particular classes of functional behavior. Some additional systems were also studied in the process of discovering these systems. Our observations for all 3-state systems are given below: 3-state systems Digraph # of free parameters classes of behavior 2 1 - (1 + b)e->-it + be- >- 2t 1 - e->-it + be->- 1 tt 3 3 3 1 - (1 + b)e->- 1 t + be->-2t 1 - e- >-it + be->- 1 tt 4 Looking at this last table, we can draw some basic conclusions about 3-state systems as a whole. We note that when we call something a unique class of functional behavior, this means that it can only be induced in a particular 3-state system and no other. Then if a class of functional behavior is common to more one system, 31 it is not unique. As for our conclusions, no 3-state system with unique functional behavior has been found, a minimum of 1 class of functional behavior and maximum of 2 classes of functional behavior always occur, and the existence of reversibility in a 3-state system leads to a failure to generate degenerative (resonance) behavior and its associated class. We hypothesize that more complex systems (e.g. reversible systems) are less likely to reduce to simple systems than simpler ones (e.g. one-way connected systems). These mathematical characterizations of functional behavior for 3-state systems will form the cornerstone of our understanding of larger reduced models and partially completes one of the goals of our research, to mathematically characterize the functional behavior of certain classes of smaller systems that we will later be studying topologically. 3.2 4-state Systems Next, we would like to mathematically characterize the functional behavior of P1i(t) for 4-state systems as we have done for 3-state systems. However, we note that we will be limited in our analyses as not all 4-state systems can be solved analytically. This limitation stems from the fact that many 4-state systems possess characteristic polynomials of order greater than 4 and thus cannot be solved analytically. In constructing the few systems that can be solved analytically, we again assume that all systems have a single absorbing final state. As for 3-state systems, we begin by enumerating all observed classes of functional behavior for a 4-state system resulting from analytically solving for their time dynamics. Note that these have been written using the minimum number of parameters based on the fact that P 41 ( t = 0) = 0: Possible classes of functional behavior (4-state systems) 1 - (1 + b + c)e->.it + be->- 2 t + ce->. 3 t 1 - (1 + b)e->.it + be->- 2 t + ce->- 2 tt 1 - e->.it + be->- 1 tt + ce->- 1 tt 2 32 Note the resonance (involving a polynomial) classes of behavior that result from degenerate eigenvalues as distinguished from the class of behavior that results from distinct eigenvalues. Next, we examine 5 simple systems under the assumptions above. Before going on to state what we observed, we give the eigenvalues Ai, ,\ 2, >.3 and the coefficients for these systems. For every system, we will provide three solutions (or state that one doesn't exist), one for the case of distinct eigenvalues with general form 1- (1 + b + c)e->- 1 t + be->-.2 t + ce->-.st, one for the case of doubly degenerate eigenvalues with general form 1-(1 +b)e- >-. 1 t+be->-. 2 t+ce->-. 2 tt, and one for the case of triply degenerate eigenvalues with general form 1 - e->-.it + be->-. 1 tt + ce->-. 1 tt 2 . Note that some systems will have no doubly or triply degenerate solution assuming positive rates. Figure 3. 7: First 4-state system. Beginning with the system depicted in Figure 3.7, we found that its distinct solution is given by the parameters: Its doubly degenerate solution is given by the parameters: Other permutations for double degeneracy yield similar expressions for b,c, as is true for the other systems we will examine. Its triply degenerate solution is given by the parameters: For the next system depicted in Figure 3.8, we found that its distinct solution is given by the parameters: A1 = r2i, A2,3 = ~(r23 + r32 + r43 ± J(r23 + r32 + r43) 2 - 4r32r43), 33 Figure 3.8: Second 4-state system. b C= ' r21r32r43 °'b,c+/3b,c+l'b,c' r21r32r43-(r21r23+r21r32+r21T43+r32T43)(r23+r32+r43±y'(r23 + r32 + r43) 2 - 4r32T43), 2 2 f3b,c = ~(r21 + r23 + r32 + T43)(r23 + r32 + r43 ± y'(r23 + r32 + r43) - 4r32r43) , o:b,c = {b,c = -~(r23 + T3z + T43 ± y'(r23 + T32 + T43) 2 - 4r32T43) 3 A doubly degenerate solution for this system has been omitted for brevity but is indeed possible. A triply degenerate solution for this system is not possible due to algebraic constraints. Figure 3.9: Third 4-state system. For the third system depicted in Figure 3.9, we found that its distinct solution is given by the parameters: Its doubly degenerate solution is given by the parameters: Its triply degenerate solution is given by the parameters: 34 Figure 3.10: Fourth 4-state system. For the fourth system depicted in Figure 3.10, we found that its distinct solution is given by the parameters: A1 = r21 + r31, A2,3 = b ~(r23 + r32 + r43 ± J(r23 + r32 + r43) 2 - 4r32r43), _ 2r21r32r43+r31r32r43 - r23r31r43-r31rl3=Fr31r43 'C - (r23+r32+r43) 2 -4r32q3 2(ab,c+.Bb,c+l'b,c ) ' ab,c = r21r32r43+r31r32r43-(r21r23+r23r31 +r21r32+r31r32+r21r43+r31T43+r32T43)(r23+ r32 + r43 ± J(r23 + r32 + r43) 2 - 4r32r43), f3b,c = ~(r21 + r23 + r31 + r32 + T43)(r23 + r32 + T43 ± J(r23 + r32 + T43) 2 rb,c = -!(r23 + T32 + T43 ± J(r23 + T32 + T43) 2 - 4r32T43) 2 4r32T43) , 3 A doubly degenerate solution for this system has been omitted for brevity but is indeed possible. A triply degenerate solution for this system is not possible due to algebraic constraints. Figure 3.11: Fifth 4-state system. Finally, for the last system depicted in Figure 3.11 , we found that its distinct solution is given by the parameters: A1-- r 21 + r 31, A2-- r 42, A3-- r 43, b - r -(rr 21 +r 42 21 31 )' c- r -(rr 31 +r ) Its doubly degenerate solution is given by the parameters: 43 21 31 . 35 Its triply degenerate solution is given by the parameters: Shown below are these 5 simple systems with their number of free parameters and classes of functional behavior. The number of classes of functional behavior for each of the digraphs was determined by studying rate-limiting conditions (e.g. Ai » A2 ) and degeneracy conditions (e.g. Ai = A2 ), and then numerically simulating that the associated classes could indeed be found: 36 4-state systems Digraph # of free parameters classes of behavior 3 1 - (1 + b + c)e->- 1t + bc>-2 t + ce->.3t 1 - (1 + b)e->- 1t + be->-2 t + ce->- 2 tt 1 - e->.it + be->- 1 tt + ce- >- 1 tt 2 4 1 - (1 + b + c)e->- 1t + be->-2 t + ce->.3t 1 - (1 + b)e->.it + be->-2 t + ce->- 2 tt , b =/=- -1 4 1 - (1 + b)e->.it + be->- 2 t + ce->- 2 tt 1 - e->.it + be->- 1 tt + ce->- 1 tt 2 5 4 From this table we see that in the group given there are a minimum of 2 classes of functional behavior and maximum of 3 classes of functional behavior, and that the existence of reversibility again leads to a failure to generate some degenerative (resonance) behavior and thus results in more limited functional behavior. Specifically, a single instance of reversibility prevents degeneracy between two of the eigenvalues (e.g. ,\ 1 =/=- >. 2 ) of the system, thus precluding the triply degenerate 37 class of behavior from occurring. However, this does not prevent degeneracy between one of these eigenvalues and another from occurring (e.g . .\ 1 = .\3 ), thus permitting the existence of the doubly degenerate class of behavior l-(l+b)e->.. 1 t+be->.. 2 t+ce->.. 2 tt with some restrictions (b =/= -l) due to rate-limiting conditions for the second and fourth digraphs in the table above. The fifth digraph is particularly interesting because though it possesses no reversibility, only doubly degenerate behavior is permitted and hence only one degenerative class of behavior occurs (in addition to the class of behavior that always occur with distinct eigenrates). It is the quintessential example that triply degenerate behavior can sometimes be prevented in systems lacking reversibility. We have found that for both 3- and 4-state systems the existence of reversibility leads to the functional form of its time dependence being more limited and, we hypothesize more generally, that the more connected a system is, the fewer classes are permissible for its functional form. Another hypothesis we propose is that as a system increases in size with only one-way connections, its ability to be reduced by means of the model reduction algorithm discussed earlier increases as a result of the greater number of opportunities for reduction. Thus this hypothesis predicts which systems will reduce well and which will not. This will be tested in future work. This concludes our mathematical characterization of 4-state systems, and reinforces some of the properties we have found for 3-state systems while adding additional depth to our understanding of smaller systems that we will later be studying topologically. 3.3 The Trapping-state System Figure 3.12: The trapping-state system. .. 38 We now take the time to study a particular 3-state system of interest. We will be examining the trapping-state system depicted in Figure 3.12. Since this system possesses reversibility, it will naturally only assume distinct eigenvalues. Recall that its distinct solution of general form 1-(1 +b)e->- 1 t +be->-2 t is given by the parameters: A1,2 = -Hr12 + r21 + r31 ± J(r12 + r21 + r31) 2 - 4r12r31), r31(r12-r21-r31+y(r12+r21+r31) 2 -4r12r31) b= 2(r12r31 +(r12+r21 +r31)(r31-r12-r21 +y1(r12+r21 +r31)2-4r12r31)+i(r12+r21 +r31-J(r12+r21 +r31) 2 -4r12r31) 2 ) We would like to explore our system's remaining degrees of freedom. To do so, we take the first step of scaling time out of our system by letting t -+ { 1 since time is common to all systems that are produced and hence we can remove the freedom due to this variable. Equivalently we may set r 31 =l. One natural way of investigating the remaining freedom in the parameters of our system is by plotting b v. fi. The corresponding equations that we will plot are then given by: ~ _ >-1 - b= l+r21+r12-y(I+r21+r12) 2-4r12 l+r21 +r12+y1(1+r21 +r12) 2-4r12' (r12-r21-l+y (ri2+r21 +1) 2-4r12) 2(ri2+(ri2+r21 +l)(l-r12-r21 +y1(r12+r21 +1) 2-4r12)+i(r12+r21 +1-J (r12+r21 +1) 2 -4r12) 2) Figure 3.13: Parameter-space of the trapping-state system. Unfortunately this equation cannot be rewritten in terms of only two variables, so we perform a 3D parametric plot with y=b, x = fi, and z=O. Since b v. fi describes an area, we will simply obtain an area in the z=O plane. I let r 12 , r 21 be 0 ... 39 (optimum values for maximum area). The plot is shown in Figure 3.13. 39 Thus we have a triangular-shaped parameter space. To more precisely characterize some of the features of the plot, we tried finding the critical points of~ v. b. Since this is a parametric function y = y(r 12 , r 21 ), x = x(r12 , r 21 ) this simply meant using the chain rule to find when £1L = ..!!JL 8r12 + ..!!JL 8r21 = 0 ax 8r12 8x 8r21 8x . Figure 3.14: Bounded area of the trapping-state system. Using numerical methods to look for solutions to this equation, no critical points were found for real, positive rates. However, rough bounds on all three sides for the area are given. These are graphically shown in Figure 3.14. The three bounding curves are b = 70~ (top left) b = 1.1(~ - 0.9) 2 + 0.06 (interior) b = 0.09~ (bottom right). The intersections of these bounding curves lie at(~, b) = (0.7959, 0.0719), (~, b) = (0.0132, 0.925), (~, b) = (0, 0). Carefully inspecting the plot suggests that b =!= 1. This has been numerically confirmed by solving the original equation b=l, yielding only negative or complex rates. Our plot more broadly tells us that given certain observed kinetics, there is only a range of possible values for those parameters that have not 40 been experimentally measured. The investigation of the parameter space of other reversible 3-state systems may suggest more insightful properties of these systems as a whole, as well as other systems. However, it is a satisfying result to have some sense of the limits of parameter space as we have found them. Thus, so far in this chapter, we have succeeded in accomplishing one of the two goals of our research, to mathematically characterize and gain intuition into the functional behavior of some of the smaller systems that will result from reducing larger systems in the remainder of our work. This understanding will prove crucial to our understanding of the signature topological features of microscopic structure that we will soon uncover and to eventually being able to make some predictions about the way in which hybridization occurs at the level of secondary structure. 3.4 Distance between Two Digraphs Having learned a good deal about how 3-, 4-, and n-state systems behave, we conclude by proposing a measure of our own for investigating the difference between two digraphs A and B. The general idea behind coming up with this measure is as follows. Degenerate behavior is what allows us to distinguish between distinct digraph topologies, so why not try to measure the distance from maximum degeneracy for a given digraph and normalize this for comparability. Then, to extend this, why not measure the difference between the normalized distances from maximum degeneracy for two digraphs and thus be able to compare how much more different one digraph is than another. In this line of thought, we offer the following measure: ~(A, B) = L (I.Ai - Ajl) - L (l.Ak - .All) ' Ai,j i \ _ /\i,j = A A.i+Ai \ _ 2 ' Ak,l = k A.1c+A.1 2 ' Ak,l B Aj is the minimizer of I.Ai - .Ajj(i #- j), .Ak is the minimizer of I.Ak - .Al I(k #- l) In words, the idea is to sum over all of the eigenvalues for each digraph the minimum absolute difference between each eigenvalue and its closest neighbor, and 41 divide by their average to normalize the value. In doing so, a small difference between eigenvalues for a large average eigenrate will be quantified as a small difference, as it should, and a small difference of comparable size between eigenvalues for a slow average eigenrate will be quantified as a more significant difference. This makes sense, since slower eigenrates should have more of a rate-limiting effect on a system and hence should be given greater weight. We hope that this measure will be useful in predicting major macroscopic changes from small microscopic changes. However, this will have to be tested in future work to see whether it bears fruit. 42 Chapter 4 The Signature Topological Features of Microscopic Structure Now that we have mathematically characterized some simple systems, we are wellprepared to address the second goal of our research, to seek out signature topological features that correspond to particular classes of macroscopic structure. The identification of these signature features will allow large, complex systems to be immediately reduced by recognition of these features within their structure and then mathematically characterized using the profiles developed in the previous chapter. And by doing so, we will lay the groundwork for understanding what small microscopic changes will lead to large macroscopic changes for our systems and thus gain greater insight into the manner in which hybridization occurs. So far, we have discovered two such signature topological features, the loop and the linear chain. These have been found by writing down the corresponding transition rate matrices A for many structurally similar digraphs and then reducing them using the model reduction scheme described earlier, and trying to generalize observations. We describe our findings here. 4.1 The Loop The loop (Figure 4.1) is a particularly interesting feature because, regardless of its size, it can be reduced to a linear 3-state system. The feature was discovered by constructing many loops of different sizes and writing down their transition-rate ma- 43 Figure 4.1: Reduction of a loop. trices and reducing them using the model reduction algorithm proposed earlier to the structure shown in the figure. The power of such a feature is that the loop effectively characterizes the opposing arrows of the reduced system, creating reverse reactions where these previously did not exist. This intuitively makes sense since the cycle that the loop represents effectively means reaching the same end state after a greater amount of time has passed (or the occupation of more states before reaching the final state), which is easily represented by adding reversible arrows to a 3-state system. Figure 4.2: Reduction of a linear chain. 4.2 The Linear Chain The linear chain (Figure 4.2) is another important feature that is often present in many complex systems. It was discovered by constructing many linear chains of different lengths and observing that their corresponding transition-rate matrices could consistently be reduced to ones characteristic of the structure shown in the figure. The idea underlying this feature is that any large linear chain can be effectively represented by a short linear chain whose rates approximate the rates of the larger system. This result is fairly intuitive and makes sense. Since such a feature will often 44 be present, the applicability of using this topological feature to simplify systems is great. While we are very excited about what we have found, we offer the caveat that these must again be confirmed once CTS form has been restored (see Figure 2.4). Such form is important for correct physical interpretation and restoration to this form has proved a formidable task for graduate student Kevin McHale (and formerly Ben Rahn) to accomplish. In the future, hopefully more signature topological features can be identified, but these two provide a solid start for exploring what small changes in microscopic structure lead to significant changes in macroscopic structure and hence understanding secondary structure kinetics theory. . 45 Chapter 5 Discussion I I · ~ I I I I Figure 5.1: Example of a large system. 5.1 Theoretical Implications Now that both some signature topological features of microscopic structure and their underlying macroscopic mathematical behavior have been established, we can begin 46 to consider the implications of our results. Consider the system shown in Figure 5. 1. The encircled regions represent portions of the system that may be simplified and this figure as a whole summarizes much of the work we have done. Since we are now able to understand both how some signature topological features are simplified and what their corresponding mathematical behavior is like, we have achieved a holistic view that will allow us to begin to think about the theory of secondary structure kinetics and to understand real systems in future work. 5.2 Future Work A number of important areas exist for future exploration. In the previous two chapters, we established the ability to predict macroscopic structure on the basis of signature topological features of microscopic structure, as well as developed an understanding of the underlying mathematical behavior of these features. We now propose a number of interesting areas for future research. Our first question is, "Can we quantify the difference between two digraphs A and Busing matrix norms, i.e. as llA- Bii?" More specifically, "Can we find a suitable metric to identify small microscopic changes that lead to significant macroscopic changes?" While we have proposed one possible measure of our own in an earlier chapter, we discuss some other possible measures here. Just as the Euclidean norm can be used to measure the distance between two vectors, we hope that a matrix norm such as the Frobenius norm, 2-norrn, 1-norm, or oo-norm will prove useful in measuring the differences between two digraphs. Tabulating these known norms alongside what our vector intuition suggests about them we have: 47 Matrix norms Measure Definition Frobenius norm llAl l} = Hypothesis L laijl = tr(A* A) crude measure .,/>:;;;;;; no intuition 2 i,j 2-norm 11 A I I2 = 1-norm llAlli = m~ L laijl = the largest absolute column sum poor measure good measure llAlloo = m~ L laijl=the largest absolute row sum J oo-norm . i i j The Frobenius norm seems the crudest norm of the bunch, resulting in a magnitudelike quantity similar to the Euclidean norm. The significance of the 2-norm is not immediately clear. The oo-norm may be interesting because it corresponds to the sum of the rates for different states whereas the 1-norm may be less so because it corresponds to the sum of rates for the same state. To test the utility of such norms for identifying small microscopic changes that lead to major macroscopic changes, we must devise pairs of digraphs whose difference leads to the loss of a signature microscopic topological feature. Then we hope to be able to use one of the above matrix norms or another measure to be able to distinguish such a change from one that leads to no loss of a signature topological feature. This work may ultimately lead to the identification of phase-transition like conditions for the systems we are studying. Thus our quest for a suitable metric is an important question that will hopefully allow us to make interesting conclusions about the way in which secondary structure kinetics proceeds and its theory as a whole. Our next question is, "How can our work be experimentally examined?" To experimentally test our work, we will need to specify what our initial conditions will be. One idea might be to engineer a strand with a particular function in mind, e.g. a nanoscale timer, since we are in such a great position to control the time dynamics, or delay, of our system. With a particular system in mind, say two specific strands of RNA, we must then generate all possible secondary structures. After we have these, we will need to construct a digraph from them and store this information in a transition-rate matrix. However, before we can do this, we will need to determine 48 the transition rates for our system. To compute the transition rates, we will need to use existing free-energy values that have been tabulated by others such as SantaLucia and Hicks [4]. Specifically, we will want to test two different measures for transition rates [10]. The first is called the Metropolis rate and is defined by where D.G = Gj - Gi· This measure makes uphill energy gradient steps more difficult, but treats downhill steps with equal probability. By contrast, the Kawasaki rate is defined by -AG rij = e 2kT • This symmetric definition accounts for the gradient in either uphill or downhill steps. The latter measure at first sight appears to be the better one because of its ability to differentiate between both uphill and downhill steps. However, both will have to be tested to be sure. Once the transition rates have been computed and the digraph reduced using the methodology outlined earlier, the next step will be to determine whether the transition rates computed hold out experimentally. This can be done using FRET (fluorescence resonant energy transfer), which allows for the measurement of nanometer scale distances between labeled sites on, e.g., RNA, and thus will permit us to study the secondary structure kinetics of RNA in realtime at the single molecule level. Such experimental work will determine whether the results of our theoretical work will be fruitful in application and hence is some of the most important work we will do in the future. Thus, it is whether our model is a sufficiently good one and whether we are able to find a way to identify what particular microscopic features play a role in significantly changing macroscopic features that will determine how well we will be able to understand secondary structure kinetics theory and apply our knowledge of it later on. 49 Bibliography [1) Gacy, A.M., Goellner, G., Juranic, N., Macura, S., and McMurray, C.T. (1995) Trinucleotide repeats that expand in human disease form hairpin structures in vitro. Cell, 81, 533-540. [2) Breaker, R.R. (1997) DNA aptamers and DNA enzymes. Curr. Opin. Chem. Biol., 1, 2631. [3) Yurke, B., Turberfield, A.J., Mills, A.P. , Jr., Simmel, F.C., and Neumann, J.L. (2000) A DNA-fuelled molecular machine made of DNA. Nature, 406, 605-608. [4) SantaLucia, J ., and Hicks, D. (2004) The Thermodynamics of DNA Structural Motifs. Annu. Rev. Biophys. Biomol. Struct., 33, 415440. [5) Rahn, B. (2001) A Balanced Truncation Primer. quant-ph/0112066. [6) Zhou, K., Doyle, J.C., and Glover, K. (1996) Robust and Optimal Control. Prentice Hall, New Jersey. [7) Dullerud, G.E. , and Paganini, F.G. (2000) A Course in Robust Control Theory: A Convex Approach. Springer, New York. [8) Gugercin, S., and Antoulas, A.C. (2004) A survey of model reduction by balanced truncation and some new results. Int. J. Control, 77, 748-766. [9) Perko, L. (2001) Differential Equations and Dynamical Systems. Springer, New York. . .. ... ' v 50 (10] Chan, H.S., and Dill, K.A. (1998) Protein folding in the landscape perspective: chevron plots and non- Arrhenius kinetics. Proteins: Structure, Function, and Genetics, 30, 2-33. ~ ~ ,. . .... ... 51 Appendix A Code for Model Reduction [t,s]=eig(A); d=diag(s); k=find(d==O); to=t (: ,k); ti=inv(t); tio=ti(k,:); po=to*tio; % construct orthogonal projector v=t; vi=ti; v(:,k)=[]; % throw out the k-th column of v vi(k,:)=[]; % throw out the k-th row of v Bt=vi*B; Ct=C*v; Dt=D; stable-111atrix=ss(At,Bt,Ct,Dt); stable-111atrix=ss(stable-111atrix,'minimal'); red_stable-111atrix=balancmr(stable-111atrix,3); ., .. 52 [a,b,c,d]=ssdata(red_stable....matrix); x=sym ( 'x ' ) ; y=c*expm(a*x)*b+C*po*B