Engineering Genetic Circuits Chris J. Myers Lecture 3: Differential Equation Analysis Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 1 / 48
Albert Einstein Yes, we have to divide up our time like that, between our politics and our equations. But to me our equations are far more important, for politics are only a matter of present concern. A mathematical equation stands forever. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 2 / 48
Introduction Next step of the engineering approach is analysis. Goal of analysis is to be able to both reproduce experimental results and make predictions in silico . Simulation provides unlimited controllability and observability. Traditional classical chemical kinetics (CCK) utilizes ordinary differential equations (ODE) to represent system dynamics. Law of mass action can translate a chemical reaction model into ODEs known as reaction rate equations . ODEs typically analyzed using numerical simulation. Qualitative analysis utilized to understand behavior as initial conditions and parameter values vary. Partial differential equations (PDE) utilized when spatial considerations are important. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 3 / 48
Overview Classical chemical kinetics Differential equation simulation Qualitative ODE analysis Spatial methods Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 4 / 48
A Classical Chemical Kinetic Model A CCK model tracks concentrations of each chemical species (i.e., number of molecules divided by volume, Ω , of cell or compartment). Assumes reactions occur in a well-stirred volume (i.e., molecules are equally distributed within the cell) and spatial effects can be neglected. Assumes reactions occur continuously and deterministically. Requires that the number of molecules of each species are large. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 5 / 48
A Classical Chemical Kinetic Model (cont) A CCK model is composed of n chemical species { S 1 , . . . , S n } and m chemical reaction channels { R 1 , . . . , R m } . Each reaction R j has the following form: k f v p 1 j S 1 + ... + v p v r 1 j S 1 + ... + v r → nj S n nj S n ← k r where v r ij is the stoichiometry for species S i as a reactant in reaction R j and v p ij is the stoichiometry for species S i as a product in reaction R j . ij and/or v p The values of v r ij are 0 when species S i does not participate as a reactant and/or product in reaction R j . Parameter k f is forward rate constant while k r is reverse rate constant . If the reaction is irreversible , then k r is 0. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 6 / 48
A Classical Chemical Kinetic Model (cont) Law of mass action states that the rate of an irreversible reaction is proportional to the product of concentrations of reactant molecules. The rate of a reversible reaction is also reduced by a value proportional to the product of the concentrations of product molecules. Formally, the reaction rate V j for reaction R j is defined as follows: n n [ S i ] v p [ S i ] v r ij − k r ∏ ∏ = V j k f ij i = 1 i = 1 where [ S i ] is the concentration of species S i . An ODE model can be constructed as follows: m d [ S i ] ∑ = v ij V j , 1 ≤ i ≤ n dt j = 1 where v ij = v p ij − v r ij (i.e., the net change in species S i due to reaction R j ). CCK model consists of one ODE for each species which is sum of the rates of change of species due to each reaction that affects the species. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 7 / 48
ODE Model Example Constant Value k d − → () RNAP 0 30 nM CI 0 . 1 M − 1 k d K d − → () CII 0 . 01 M − 1 K o 1 K o 1 P RE + RNAP ← → S 1 0 . 69422 M − 1 K o 2 K o 2 P R + RNAP ← → S 2 0 . 2165 M − nc K r k b 0 . 00161 M − ( na + 1 ) − → S 1 + np CI S 1 K a 0 . 014 sec − 1 k o k o − → S 2 + np CII S 2 0 . 00004 sec − 1 k b K d ← → 2CI CI 2 0 . 015 sec − 1 k a K r P R + nc CI 2 ← → S 3 0 . 0075 sec − 1 k d K a P RE + na CII + RNAP ← → S 4 nc 1 k a − → S 4 + np CI na 1 S 4 np 10 Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 8 / 48
ODE Model Example d [ CI ] np k b [ S 1 ]+ np k a [ S 4 ] − 2 ( K d [ CI ] 2 − [ CI 2 ]) − k d [ CI ] = dt d [ CI 2 ] K d [ CI ] 2 − [ CI 2 ] − nc ( K r [ P R ][ CI 2 ] nc − [ S 3 ]) = dt d [ CII ] np k o [ S 2 ] − na ( K a [ P RE ][ RNAP ][ CII ] na − [ S 4 ]) − k d [ CII ] = dt d [ P R ] [ S 2 ] − K o 2 [ P R ][ RNAP ]+[ S 3 ] − K r [ P R ][ CI 2 ] nc = dt d [ P RE ] [ S 1 ] − K o 1 [ P RE ][ RNAP ]+[ S 4 ] − K a [ P RE ][ RNAP ][ CII ] na = dt d [ RNAP ] = [ S 1 ] − K o 1 [ P RE ][ RNAP ]+[ S 2 ] − K o 2 [ P R ][ RNAP ]+ dt [ S 4 ] − K a [ P RE ][ RNAP ][ CII ] na d [ S 1 ] = K o 1 [ P RE ][ RNAP ] − [ S 1 ] dt d [ S 2 ] = K o 2 [ P R ][ RNAP ] − [ S 2 ] dt d [ S 3 ] K r [ P R ][ CI 2 ] nc − [ S 3 ] = dt d [ S 4 ] K a [ P RE ][ RNAP ][ CII ] na − [ S 4 ] = dt Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 9 / 48
Differential Equation Simulation The differential equations for a set of reaction rate equations are: dx i = f i ( x ) , where 1 ≤ i ≤ n dt where x = [ x 1 ,..., x n ] ≥ 0 is vector of species concentrations. Solving this ODE model analytically is very difficult, if not impossible. Numerical simulation can approximate time evolution of X ( t ) assuming X ( t 0 ) = x 0 ( initial value problem ). Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 10 / 48
Euler’s Method Simplest approach to solve the initial value problem is Euler’s method . Initial instantaneous rate of change for each S i at time t 0 : dX i ( t 0 ) = f i ( x 0 ) , where 1 ≤ i ≤ n . dt If the rate of change, f i ( X ( t )) , remains constant for all t ≥ t 0 , then X i ( t ) = x 0 i + f i ( x 0 )( t − t 0 ) . Not true in general, but may be reasonable to assume value remains close to f i ( x 0 ) for some small time step ∆ t ( step size ). Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 11 / 48
Forward Euler Method With this assumption, X i ( t 1 ) ≈ X i ( t 0 )+ f i ( X ( t 0 )) ∆ t In general, for any t j = t 0 + n ∆ t where j = 0 , 1 , 2 , 3 ,... : X i ( t j + 1 ) ≈ X i ( t j )+ f i ( X ( t j )) ∆ t This algorithm is known as the forward Euler Method , and it is an example of an explicit ODE simulation method. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 12 / 48
Forward Euler Method X i ( t ) t 0 t 1 t 2 Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 13 / 48
Backward Euler Method Backward Euler method is an implicit ODE simulation method. New rate of change cannot be determined directly from current state, but it must be found implicitly by an equation that must be solved. Backward Euler method is defined by follows: X i ( t j + 1 ) ≈ X i ( t j )+ f i ( X ( t j + 1 )) ∆ t . New value of X ( t j + 1 ) is determined as the rate at the point that would have taken you there in a ∆ t step. Since X ( t j + 1 ) is not yet known, this equation must be solved using a root finding technique such as the Newton-Raphson method. Implicit methods are more complicated, but often more stable for stiff equations (i.e., those that require a very small time step). Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 14 / 48
Backward Euler Method X i ( t ) t 0 t 1 t 2 Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 15 / 48
Drawback of Euler’s Method Using the fundamental theorem of calculus, exact solution for each species, S i , must satisfy: Z t j + 1 X i ( t j + 1 ) = X i ( t j )+ f i ( X ( t )) dt t j Drawback of Euler methods is they approximate this integral by assuming that f i ( X ( t )) is constant throughout the entire ∆ t interval from t j to t j + 1 . Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 16 / 48
Midpoint Method (Second-Order Runge-Kutta) Approximates rate of change in ∆ t interval using rate at the midpoint: � � f i ( X ( t j )+ 1 X i ( t j + 1 ) ≈ X i ( t j )+ 2 f ( X ( t j )) ∆ t ) ∆ t X i ( t ) t 0 t 1 t 2 Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 17 / 48
Fourth-Order Runge-Kutta More points can be combined to further improve accuracy: α 1 = f ( X ( t j )) � X ( t j )+ ∆ t � α 2 = 2 α 1 f � X ( t j )+ ∆ t � α 3 2 α 2 = f α 4 = f ( X ( t j )+ ∆ t α 3 ) X ( t j )+ ∆ t X ( t j + 1 ) = 6 [ α 1 + 2 α 2 + 2 α 3 + α 4 ] Implicit Runge-Kutta methods can also be used for stiff equations. Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 18 / 48
Adaptive Stepsize Control To obtain good results, should modify ∆ t during the simulation. Simulation should slow down when rates are changing rapidly and speed up when the rates are changing slowly. With the step doubling approach, X ( t j + 1 ) , is found in one ∆ t step and ′ ( t j + 1 ) is found by taking two ∆ t / 2 steps. X Estimate of the error for each species: ′ E i = | X i ( t j + 1 ) − X i ( t j + 1 ) | Goal of adaptive stepsize control is to achieve a desired error level: D i = abs + rel ·| X i ( t j + 1 ) | where abs is absolute error level and rel is relative error level . Chris J. Myers (Lecture 3: ODE Analysis) Engineering Genetic Circuits 19 / 48
Recommend
More recommend
