Control Theory for Biologists
Contents
1. 1 2 Models and Systems 2i s 22 4 4 L3 Model Variables and Parameters 1 4 Classification of Models 2 5 2 4 4 ae4 1 5 Electtical M dels s s a coe e ee OR 16 Block Diagrams osos s s ossos eion m a paas 1 7 Digital Computer Models 2 224824 1 8 Dimensions and Units 1 9 Approximations oa o4eo44 oae bas amp 110 Be ebhaviof s se ests me GADD ABD OR S 1 11 Software for Solving ODEs 4 44 14 Exercises 00 20 00 eee eee rarr References History aon wn _ 14 17 26 28 30 32 34 35 41 47 49 iii iv CONTENTS Chapter 1 Introduction to Modeling 1 1 Build a Simulation Model Water Tank Model Figure 1 1 shows two water tanks The first tank is fed with water at arate Q m s This tank drains into a second tank at a rate Q2 which in turn drains to waste at a rate Q3 The second tank has an additional feed of water flowing in at a rate Q4 The height of the water level in each tank is given by h and hz respectively Each tank has a cross sectional area A The rate of change in the volume of water in a given tank is the rate at which the water enters minus the rate at which it leaves For example for the first tank we have dV a An If we want the equation in terms of the rate of change of height then we needs to recall that V Ah that is dV dh 2 aA dt dt 2 CHAPTER 1 INTRODUCTION TO MODELING so that dh _ Qi Q2. In this case the model is simple enough that we can derive the ana lytical solution for this model Assuming that at time zero both S2 and S3 are at zero concentration then it can be shown that S1 t S1 0 e k SO S10 e e koe pe so s0 1 2 In the majority of cases it is not possible to derive such an analytical solution and we must therefore fall back on using a computational model To do this we must find a suitable software application that can solve differential equations The plot shown in Figure 1 16 was obtained from Jarnac www sbw app org but similar plots could be obtained using Matlab 30 CHAPTER 1 INTRODUCTION TO MODELING Simulation Result S52 Si 83 T T T Species Concentrations Figure 1 16 Computer Simulation of two Consecutive Chemical Reactions 0 10 ky 1 5 k2 0 4 Assignment 4 Obtain the same plot shown in Figure 1 16 using Matlab or Python 1 8 Dimensions and Units Variables and parameters that go into a model will be expressed in some standard of measurement In science the recognized standard for units are the SI units These include units such as the meter for length kilogram for mass second for time Joules for energy kelvin for temperature and the mole for amount The mole is of particular importance because it is a means to measure the number 1 8 DIMENSIONS AND UNITS 31 of particles of substance irres
3. A v4 v_ where A is the so called gain of the amplifier In real op amps the gain is usually of the order of at least 10 For an op amp operating with a supply voltage of 15 volts the difference in voltage between the two input terminals cannot be more than 120 uV in order to avoid saturation of the amplifier This means that an op amp on its own is of little practical value because the slightest difference in voltage at the input terminals will cause the op amp to saturate When applied to real problems op amps are almost always used with some kind of feedback circuit Output Voltage Saturation Figure 1 10 Basic Op Amp Response Although the basic op amp is essentially a very high gain amplifier it has certain other properties which make it extremely useful These properties are listed below 22 CHAPTER 1 INTRODUCTION TO MODELING 1 The resistance or impedance into the input terminals is ex tremely high 2 The resistance or impedance of the output terminal is ex tremely low The high impedance on the input terminals means than an op amp will impose hardly any load on the circuit it is connected to The input resistance will often be of the order of a least a mega ohm often more A low impedance on the output terminal of the order of 70 ohms or less means that attaching loads to the op amp will have hardly any effect on the output voltage As a result of these properties op amps are ideas c
4. The complete list of rates is given below v2 koDp v3 k3D v4 kD vs ks D ve ke Db where Dp is the concentration of drug in the blood stream 1 11 SOFTWARE FOR SOLVING ODES 45 D the concentration of drug in the liver and Da the concen tration of drug at the site of action Using the conservation of mass three differential equations can be written down that describe the rate of change of drug in each of the three com partments Previous pharmokinetic studies of the drug indicate that the rate constants for the entry and exit of drug to and from the various compartments are given by k2 6 12 k3 0 2 k4 0 45 ks 1 k 5 All units in seconds Using a simulation model answer the following question Assume that the nurse can inject 1 mM of drug per second into the patient s blood stream That is in one second the concentration of drug in the blood stream increases by 1 mM A little unrealistic perhaps This is another way of saying that v 1 mM sec The nurse can start and stop the injection very quickly such that the profile one observes in the rate at which drug enters the blood looks like a pulse as shown in Figure 1 23 Rate of Infusion 1 6 Time Secs Figure 1 23 Drug Infusion Profile Assume that in order for the drug to be effective the concen tration at the site of action must reach at least 2 5 mM If the 46 CHAPTER 1 INTRODUCTION TO MODELING concentration o
5. cannot take on all values within a given numeric range For example the number of aeroplanes in the sky at any one time is a discrete number In statistics this is generalized further to a finite set of states such as true false or combinations in a die throw Continuous variables can assume all values within a given nu meric range For convenience we will often represent a mea surement as a continuous variable For example we will often use a continues variable to represent the concentration of a so lute as it is unwieldy to refer to the concentration of a solute as 5 724 871 927 315 193 634 656 molecules per liter The reason for this is that each step in the simulation is determined by one or more random processes To give an example modeling lion predation on the Serengeti could be modeled as a stochastic process It is not guaranteed that a Lion will catch its prey every time instead there is a probability it will succeed To model this process the computer simulation would throw a die to determine whether the Lion had succeeded or not Repeatedly running such a simulation again would naturally give a slightly different outcome because the die throws would be different in each run 1 3 MODEL VARIABLES AND PARAMETERS 13 A deterministic model is one where a given input will always produce the same output For example in the equation y x setting x to 2 will always yield the output 4 A stochastic model is one where th
6. include more of the environment It will often be the case that the interface between the environment and the system will not be perfect so that there will be some effect that the system has on the environment So long as this effect is small we can assume that the environment is not affected by the system Replacing complex subsystems with lumped or aggregate laws Lumping subsystems is a commonly used technique in simplifying cellular models The most important of these is the use of aggre gate rate laws such as Michaelis Menten or Hill like equations to 34 CHAPTER 1 INTRODUCTION TO MODELING model cooperativity Sometimes entire sequences of reactions can be replaced with a single rate law Assuming simple linear cause effect relationships In some cases it is possible to assume a linear cause effect between an enzyme reaction rate and the substrate concentration this is especially true when the substrate concentration is below the Km of the enzyme Linear approximations make it much easier to understand a model Physical characteristics do not change with time A modeler will often assume that the physical characteristics of a system do not change for example the volume of a cell the values of the rate con stants or the temperature of the system In many cases such approx imations are perfectly reasonable Neglecting noise and uncertainty Most models make two impor tant approximations The first is that noise in the system i
7. the rate law the rate constants in simple first order kinet ics are expressed in per unit time t while in second order reac tions the rate constant is expressed per concentration per unit time mol 71 32 CHAPTER 1 INTRODUCTION TO MODELING In dimensional analysis units on the left and right hand sides of ex pressions must have the same units or dimensions There are cer tain rules for combining units when checking consistency in units Only like units can be added or subtracted thus the expression S k cannot be summed because the units of S are likely to be mol and the units for k t Even something as innocent looking as 1 S can be troublesome because S has units of concentration but the constant value 1 is unit less Quantities with different units can be multiplied or divided with the units for the overall expression computed using the laws of exponents and treating the unit symbols as variables Class Question 7 Determine the overall units for the expression k S Km where the units for each variable are k t 1 S mol and K mol 17t In exponentials such as e the exponent term must be dimension less or at least the expression should resolve to dimensionless thus e t is permissible but e is not if for example k is a first order rate constant Trigonometric functions will always resolve to dimen sionless quantities because the argument will be an angle which can always be expre
8. these passive circuits suffer from loading issues and are not particular practical as compu tational circuits Avoiding the Loading Issue The way to avoid the loading issues described in the last section is use some kind of buffer circuit such as an Op Amp Op Amps or nttp www youtube com wat ch v JbpR5nGudGg 3http www youtube com watch v qtrYd0uJzyA 20 CHAPTER 1 INTRODUCTION TO MODELING operational amplifiers are the workhorse of many analog circuits they can be purchased very cheaply with varying characteristics to suit a myriad of different applications Figure 1 8 Common 741 Op Amp Image from Jameco Elctron ics cost 0 25 each In their basic form an op amp has five terminals Two of these ter minals are used to supply power to the circuit and are not normally shown when drawing op amp circuits The remaining three termi nals represent two input and one output terminal The two input terminals are designated the inverting and non inverting terminals respectively In diagrams the two input terminals are often indicated with a plus and a minus sign Figure 1 9 Vs Vout Vs Figure 1 9 Symbol for an Op Amp Vs and Vs are the positive and negative supply voltages V and V are the inverting and non inverting inputs 1 5 ELECTRICAL MODELS 21 Basic Properties of an Op Amp An Op Amp amplifies the difference in voltage between the positive and negative input terminals that is Vo
9. 2 dt A Using Torrielli s Law we know that the rate of water flowing out of a given tank 7 is equal to Qi K Vvh Where K is a constant related to the resistance of the output pipe Therefore for the first tank we have dh Q Kvh dt A Figure 1 1 Water Tank Model Given this information answer the following questions a Plot the rate of outflow Q2 as a function of the height of water h at a given resistance K4 b Assuming that Q and Q are fixed and we start with both tanks empty what do you expect to happen over time as water flows in c Write out the differential equations ODEs that describe the rate of change in the tank water levels h and h2 1 2 MODELS AND SYSTEMS 3 d Build a computer model of the tank system assign suitable values to the parameters in the model and run the simulation to plot the height of water in the tanks over time Assume both tanks are empty at time zero e Investigate the effect of increasing and decreasing the resistance parameters K and K3 on the model 1 2 Models and Systems In the previous section we built a quantitative model of a water tank system but what exactly is a model A model according to the Oxford English Dictionary is A simplified or idealized description or conception of a particular system situation or process often in math ematical terms that is put forward as a basis for theo retical or empirical understanding o
10. A KVL Around a cycle gt vk 0 1 5 ELECTRICAL MODELS 17 1 5 Electrical Models We are now going to briefly look at electrical systems and how they can be used to model physical systems Class Question 4 Given the following resistor network determine the voltage V Figure 1 3 Resistor Network When you ve derived the relation for vo assume R2 gt gt R and then infer what calculation could be done with the circuit Figure 1 4 Block that computes the average of inputs V to Vy Class Question 5 Show that the resistive circuit in Figure 1 6 can be used to compute a voltage times a constant 18 CHAPTER 1 INTRODUCTION TO MODELING Figure 1 6 Circuit to multiply an input voltage by a constant Class Question 6 What is the problem with all these circuits that could make them potentially poor at computing There are many other kinds of circuits one can build using passive components Two of the most common are the integrator and differ entiator circuits Figure 1 7 Integrators carry out the mathematical operation J vdt 1 5 ELECTRICAL MODELS 19 while differentiators carry out the operation dv dt a R b C Sv Vo vu Vo T r Integrator Differentiator Figure 1 7 Integrator and differentiator circuits Assignment 1 Watch the YouTube videos RC Integrator Circuit and RC Differentiator Circuit Figure 1 7 As with the summer and multiplier circuits
11. Control Theory for Biologists Herbert M Sauro University of Washington Seattle WA Ambrosius Publishing Copyright 2111 2013 Herbert M Sauro All rights reserved First Edition version 0 2 Published by Ambrosius Publishing and Future Skill Software www analogmachine org Typeset using BIEX 2 TikZ PGFPlots WinEdt and 11pt Math Time Professional 2 Fonts Limit of Liability Disclaimer of Warranty While the author has used his best efforts in preparing this book he makes no representations or war ranties with respect to the accuracy or completeness of the contents of this book and specifically disclaim any implied warranties of merchantability or fitness for a particular purpose The advice and strategies contained herein may not be suitable for your situation Neither the author nor pub lisher shall be liable for any loss of profit or any other commercial dam ages including but not limited to special incidental consequential or other damages No part of this book may be reproduced by any means without written permission of the author ISBN 10 x xxxx xxxx x ebook ISBN 13 xxx x xxxxxxx x x ebook ISBN 10 x xxxxxxx x x paperback ISBN 13 xxx x xxxxxxx Xx x paperback Printed in the United States of America Mosaic image modified from Daniel Steger s Tikz image http www texample net tikz examples mosaic from pompeii Contents 1 Introduction to Modeling tell Build a Simulation Model 4 4 4 4x
12. Lumped and Distributed Parameter Models Many complex models can be approximated with a single number For example we often describe a resistor using a single value its re sistance In reality the resistor has a length a diameter and a chem ical composition The resistance is a function of all these properties that make up the resistor We could model the resistor by slicing up the resistor in to many small compartments and compute the resis tance as a systemic property In the former case we have what is called a lumped parameter model in the second case a distributed parameter model Linear and Nonlinear Systems We will discuss linear and no linear systems in much more detail in Chapter but we can say at this stage that a linear system is one where the output changes in proportion to the inputs That is doubling the input will result in a doubling of the output 16 CHAPTER 1 INTRODUCTION TO MODELING Box 1 0 Review of basic electricity v Voltage i Current R Resistance C Capacitance Q Charge Ohms Law v R 1 Capacitor Law Q Cv or v ae In addition the charge on a capacitor is given by q iat combining the two previous relations and differentiating yields dv See dt Resistors in series parallel In series R gt R In parallel 1 R 1 R Capacitors in series parallel In series 1 C we E In parallel C ya C Kirchhoff s Laws KCL and KVL KCL At a node Dain Ss eye
13. a A projectile fired from a cannon Figure 1 20 Figure 1 20 Projectile fired from a cannon b An RC circuit acting as a low pass filter Figure 1 21 c An ecosystem comprising of a population of rabbits and foxes where the model is base on the Lotka Volterra system of equations d A chemical system based on the Brusselator model 2 Build and simulate the following models 42 CHAPTER 1 INTRODUCTION TO MODELING Ry NNT Vo C Figure 1 21 RC Circuit Predator Prey Models a Build a computer model of the predator prey system de scribed by the Lotka Volterra system of equations dP a Prey p Prey Predator dPredat Prey Predator y Predator where is the growth rate of the prey 6 is the rate of pre dation y is the death rate of the predator and is the rate of growth in the predator population b Explore how the parameters affect the evolution of the pop ulations Biochemical Models Build a computer model of the following simple metabolic pathway vi v2 v2 Xo Si gt So gt X 1 11 SOFTWARE FOR SOLVING ODES 43 Assume that the metabolites Xo and X are fixed Assume also that v and v2 are governed by the rate laws vy Ey ky Xo k2S V2 Ex k3S k4S2 v3 E3 ks So where is the concentration of the i enzyme and k are the kinetic rate constants a Write out the variables and parameters of the system b Wri
14. block to consider it the integral block The input to this block will be the derivative dS dt and the output will be S4 How ever S appear as a negative term in the differential equation k S4 therefore we will also invert S 2 at the same time Once we have S we only need to add a multiply by a constant block to multi ply S by the constant k to obtain the value of dS dt This final operation closes the loop Figure 1 15 illustrates the complete block diagram From the block diagram we could construct this device out of suit ably sired op amps such that the electrical circuit would solve the differential equation Before the advent of the digital computer such analog computers were common place and were use to solve a va riety of problems in science and engineering One of the advantages of an analog computer is that they are very fast and can be oper ated in a life mode such that it is possible to interact with the model in real time and the analog computer will automatically adjust the solution However with the advent of fast and particularly cheap digital com puters the need for analog computer has declined although there is increasing interest to use some of the ideas in a new field called 1 6 BLOCK DIAGRAMS 27 y Summation yo JP o Y2 Y3 Y3 Yo Integration J gt Y Mult by Constant y gt ky Multiply Variables 41 Yr Y2 Yo Figure 1 14 Common Block Symbols Neuromorhpic electron
15. d 4 January 2013 8 Wikipedia 2013d http en wikipedia org wiki Difference_engine Online accessed 4 January 2013 9 Wikipedia 2013e http en wikipedia org wiki Differential_analyser Online accessed 4 January 2013 10 Wikipedia 2013f http en wikipedia org wiki Emergence Online accessed 4 January 2013 11 Wikipedia 2013g http en wikipedia org wiki Fractal On line accessed 4 January 2013 12 Wikipedia 2013h _ http en wikipedia org wiki MONIAC _Computer Online accessed 4 January 2013 13 Wikipedia 20131 http en wikipedia org wiki Neural_networks Online accessed 4 January 2013 14 Wikipedia 2013j http en wikipedia org wiki Rangekeeper Online accessed 4 January 2013 15 Wikipedia 2013k _ http en wikipedia org wiki Slide_rule Online accessed 4 January 2013 16 Wikipedia 20131 http en wikipedia org wiki SPICE On line accessed 7 January 2013 17 Wikipedia 2013m http en wikipedia org wiki Tide _predicting machine Online accessed 4 January 2013 18 Wikipedia 2013n http en wikipedia org wiki V 2a Online accessed 4 January 2013 History 1 VERSION 0 1 Date 2011 04 2 Author s Herbert M Sauro Title Introduction to Control Theory for Biologists Modification s First Edition Release 49 50 BIBLIOGRAPHY
16. e observer or under the control of a forcing function There is a third set of quantities which are often called the system parameters These include the various kinetic constants and enzyme activity factors whose values are determined by a com bination of thermodynamic physical or genotypic properties Some system parameters will appear in the reaction rate laws as kinetic constants others will be related to physical characteristics such as the volume of the system assuming the volume is constant during the study period As with boundary variables the system param eters are in principle under the control of the experimenter and are not a function of the model itself One can imagine for example 1 3 MODEL VARIABLES AND PARAMETERS 11 changing kinetic constants via site directed mutagenesis or chang ing enzyme activities by altering promoter efficiencies The actual choice of system parameters in any particular model will depend on physical limitations and the question being considered by the re searcher Class Question 2 For the water tank model described in section 1 1 write out the vari ables and parameters of the system Identify the boundary parame ters Mathematical Descriptions of Models Table 1 1 lists the many different ways in which a model can be represented and simulated However one thing they all have in com mon is that the models are expressed first in mathematical form The form of this expression determin
17. e processes described by the model include a random element This means that repeated runs of a model will yields slight different outcomes We can now therefore classify a model as a combination of the above attributes The water tank model uses a deterministic continuous approach The lion model might use a discrete and stochastic ap proach In this book we will be primarily concerned with determin istic and continuous models Table 1 2 shows the four combinations and example where each combination might be appropriately used Type Example Continuous Deterministic Projectile motion Continuous Stochastic Brownian motion Discrete Deterministic Large population dynamics Discrete Stochastic Small population dynamics Table 1 2 Examples of different kinds of model Class Question 3 How you would describe each of the following systems in terms of how you might model them A game of chess Two teams playing American Football A forest fire Gene regulation in a bacterium The beating of a healthy heart ofc a ee 14 CHAPTER 1 INTRODUCTION TO MODELING A tornado Growth of a tumor Upper atmosphere chemistry A digital computer 10 Electron flow through a single transistor 11 A population of cells from which emerges a cancerous cell 12 Metabolism 13 The swimming of a single E coli through water Oo OND 1 4 Classification of Models In addition to classifying models as discrete continuous and detem
18. es how the model will be con structed and how it will be solved or simulated In particular deci sions must be made about how the variables and parameters in the system are best described For example the model variables can be described either using discrete or continuous variables For exam ple the change in the level of water in a tank is more reasonably described using a continuous variable such as the height On the other hand it might be more realistic to describe the dynamics of lion predation on the Serengeti using a discrete model where indi vidual lions can be represented It might not make much sense to refer to 34 67 lions in a model The choice of whether to use a dis crete or continuous description depends entirely on the system being studied and the questions posed Another important categorization is whether the model should be represented in a deterministic or stochastic form A deterministic model is one where if we repeated the simulation using the same starting conditions we would get exactly the same result again That 12 CHAPTER 1 INTRODUCTION TO MODELING is the future state of the model is completely determined by its initial starting point The model of the water tanks filling up is an example of a deterministic model In contrast each time we ran a computer simulation of a stochastic model we would get a slightly different outcome even though the starting conditions are the same A discrete variable is one that
19. f drug in the syringe is 1 mM use the com puter model to roughly estimate the minimum time the injec tion should last in order for the drug to be effective Bibliography 1 Boahen K 2005 Scientific American 292 5 56 63 2 Hucka M A Finney H M Sauro H Bolouri J C Doyle H Kitano A P Arkin B J Bornstein D Bray A Cornish Bowden A A Cuellar S Dronov E D Gilles M Ginkel V Gor I I Goryanin W J Hedley T C Hodgman J H Hofmeyr P J Hunter N S Juty J L Kasberger A Krem ling U Kummer N Le NovAASAAtre L M Loew D Lucio P Mendes E D Mjolsness Y Nakayama M R Nelson P F Nielsen T Sakurada J C Schaff B E Shapiro T S Shimizu H D Spence J Stelling K Takahashi M Tomita J Wagner and J Wang 2003 Bioinformatics 19 524 531 3 Sauro H M 2000 Jn Hofmeyr J H S J M Rohwer and J L Snoep ed Animating the Cellular Map Proceedings of the 9th International Meeting on BioThermoKinetics Stellen bosch University Press 4 Trafton A 2012 http www mit edu newsoffice 201 1 brain chip 1115 html 5 Wikipedia 201 3a http en wikipedia org wiki Analog_computer Online accessed 4 January 2013 6 Wikipedia 2013b http en wikipedia org wiki Antiky thera_mechanism Online accessed 4 January 2013 47 48 BIBLIOGRAPHY 7 Wikipedia 2013c http en wikipedia org wiki Cellular_automaton Online accesse
20. ffect of changes in gene expression The separation between the internal and external variables depends ultimately on practical considerations and the particular questions that the experimenter and modeler wish to answer However once the choice is made the separation is strictly adhered to during the course of a study This means for example that the environment sur 10 CHAPTER 1 INTRODUCTION TO MODELING rounding the physical system will by definition be unaffected by the behavior of the system If for some reason parts of the environ ment do change then these parts should be considered part of the system Forcing Functions As described earlier it is common to make sure that the surround ings do not change during the duration of the study For example we might make sure that the pH remains constant by using a buffer so lution The key point is that the experimenter has complete control over the experiment In some cases it is useful for an experimenter to change the surrounding conditions in a controlled fashion For example he she might slowly increase the concentration of an ad ministered drug or make a step change in a variable such as the con centration of an enzyme In systems theory such controlled changes are often called forcing functions Model Parameters As described previously the state variables are quantities that evolve in time and the boundary variables are concentrations or voltages that are clamped by th
21. gend S1 S2 Figure 1 18 A Matlab script for solving a simple set of differen tial equations and plotting the results Solving the Lorenz Attractor In this section we will illustrate the use of Matlab and Python to solve the famous chaotic model called the Lorenz system The 38 CHAPTER 1 INTRODUCTION TO MODELING NDSolve 3 4 2 3 xi t x1 t 2 3 x1 t 0 4 x2 t x2 t xif0 0 x2 0 0 x1 x2 t 0 5 Plot Evaluate x1 t x2 t s t 0 5 PlotRange gt A11 Figure 1 19 Mathematica script for solving a simple set of dif ferential equations Lorenz system is a three variable system of differential equations given by the following Oyen dt D px y yz det P LFY d Tafa The first example illustrates how Matlab can be used to solve the differential equations To solve differential equations we must first declare a function that computes the right hand sides of the differen tial equations and return the derivatives to the caller This is shown below function xprime lorenz t x Computes the derivatives of the differential equations sig 10 beta 8 3 rho 28 xprime sig x 1 sig x 2 rho x 1 x 2 x 1 x 3 betatx 3 x 1 x 2 1 11 SOFTWARE FOR SOLVING ODES 39 To solve the Lorenz system one of the Matlab provided ODE solver routines can be used In the script below ode45 is used which imple ments an adaptive step size Runge Ku
22. gs Any mass in the system is isolated from the surroundings An example of a close system is the bomb calorimeter Isolated System Closed System u S K Energy gt awT 2 anh Open System vry Energy se S Mass gt T Q TL Figure 1 2 Open and Closed Systems Class Question 1 For each of the following systems decide whether the system is iso lated closed or open Comment on the nature of the surroundings i A system represented by a mechanical clock slowly winds down in a room controlled by a thermostat ii A car engine running idle in the open air iii A bacterial culture is grown in batch and kept in a sealed and insulated chamber What Makes a Good Model What makes a good model There are a range of properties that a good model should have but probably the most important are accu racy predictability and falsifiablity 6 CHAPTER 1 INTRODUCTION TO MODELING gt A model is considered accurate if the model is able to describe current experimental observations that is a model should be able to reproduce the current state of knowledge gt A predictive model should be able to generate insight and or pre dictions that are beyond current knowledge Without this ability a model is considerably less useful some would even suggest useless gt Finally a model should be falsifiable By this we mean that a model cannot be proved be true only disproved The only discipline where stateme
23. he past to construct and simulate models Table 1 1 Different ways to construct and solve physical models Electrical Circuits General purpose analog computer 5 8 CHAPTER 1 INTRODUCTION TO MODELING WWII V2 guidance system 18 Neuromorhpic electronics 1 4 Mechanical and Fluid Slide rule 15 Curta 8 Tide predicting machine 17 Computing projectile trajectories 14 Differential analyzer solves ODEs 9 Antikythera mechanism planetary motion 6 Water tanks MONIAC economic model 12 Purely Mathematical Algebraic Equations Linear differential equations Linear difference equations Partial differential equations Probabilistic models Statistical models Digital Computer Solving ODEs and PDEs Agent based models multicellular systems Cellular automata 7 Emergent systems Ant models 10 Fractal models 11 Neural networks 13 1 3 Model Variables and Parameters Model Variables One of the first decisions to make when developing a model is to identify the external the causes and internal variables effects The internal variables of the system are often termed the state vari 1 3 MODEL VARIABLES AND PARAMETERS 9 ables or dynamic variables and will change in time as the system evolves The values of the state variables are determined by the model or system In an experiment involving living cells the state variables will often be concentrations of molecules and ions volt ages across membranes
24. heory of common descent is falsifiable If an entire population of fos silized rabbits were found in a precambrian formation the theory would be in serious trouble 1 2 MODELS AND SYSTEMS 7 by both the popular media and lay population There are other attributes of a model that are desirable but not es sential these include parsimonious and selective A parsimonious model is a model that is as simple as possible but no simpler Oc cam s infamous razor states that Entities should not be multiplied beyond necessity and argues that given competing and equally good models the simplest is preferred Finally since no model can rep resent everything in a given problem a model must be selective and represent those things most relevant to the task at hand Steps in Building Model To summarize we can break down the approach to building a model into at least four stages Define the system boundaries Define the simplifying assumptions Invoke physical laws to describe the system processes Test validate the model against experimental data Different Ways to Represent Physical Models The tank model described in section 1 1 was built using a set of ordi nary differential equations ODEs and solutions to these equations were obtained using a digital computer There are however many other ways to build and solve a physical model Table 1 1 lists some of the more common and interesting approaches that people have used in t
25. hey work will be given in Chapter Figure 1 18 illustrate a simple Matlab script for solving a set of two ordinary differential equations Fig ure 1 19 shows the same simulation but using Mathematica instead Both Matlab and Mathematica are however expensive to purchase therefore many researchers have either written their own solvers or use publicly available software libraries There are free Matlab like tools such as Octave SciLab Euler and Yorick that can be successfully used There is is also consider able support for solving ODEs from the popular scripting language Python In addition to these there are specialized tools for mod eling specific types of physical systems For example Jarnac has 4nttp www gnu org software octave Shttp www scilab org Shttp euler sourceforge net http yorick sourceforge net Shttp www python org www sbw app org 1 11 SOFTWARE FOR SOLVING ODES 37 specific support for describing chemical or biochemical networks Solve a systems of ODE s Use ODE45 Model gt S1 gt S2 gt using simple mass action kinetics function dxdt odes t x dxdt 1 3 4 2 3 x 1 dxdt 2 2 3 x 1 0 4 x 2 return Set an error tolerance options odeset RelTol 1e 6 Initial conditions Xo 0 0 Timespan tspan 0 5 Call the solver t X ode45 odes tspan Xo options Plot the results plot t x 1 t x 2 le
26. ical problem in a way that will be more familiar to the user From such input the software will derive the differential equations automati cally This can help reduce errors in the final equations Examples of such tools include SPICE 16 for modeling electrical circuits and Jarnac 3 for modeling biochemical networks In biochemical modeling there are also standard file formats for storing models The most well know of these is the Systems Biology Markup Language SBML 2 This means that applications that support SBML can access the large number of models available in model repositories such as Biomodels Many of the example models in this book will be described using Jarnac scripts which can be easily converted to SBML or Matlab Jarnac is a windows application that offer a wide range of analysis methods for modeling biochemical networks A simple example is given here Define a simple linear chain of reactions p defn cell Xo gt S1 k1 Xo k2 S1 S1 gt S2 k3 S1 k4 S1 S2 gt X1 k5 S2 2 end Initialize parameters p Xo 10 p ki 0 6 p k2 0 4 1 11 SOFTWARE FOR SOLVING ODES 41 p k3 p k4 p k4 0 1 3 4 1 2 Carry out simulation and plot results Three arguments timeStart timeEnd number0fPoints m p sim eval 0 10 100 graph m A user manual can be found at www sbw app org Exercises 1 Identify the state variables and parameters in the following systems
27. ics 1 4 where the aim is to recreate brain like computational systems from electrical analog circuitry Assignment 3 Draw a block diagram of the following system of two differential equa tions dS ol h8 dt 191 dS gt e dt 1921 2 28 CHAPTER 1 INTRODUCTION TO MODELING S dS dt FA k Figure 1 15 Analog System to Solve das k S Si t S1 0 eo ky kit kot S2 t Ba e e 1 7 Digital Computer Models We now turn to the use of the programmable digital computer for modeling physical systems Modern digital computers are programm able machines that can manipulate binary numbers As such they are very adept at handling problems involving numerical process ing There are a great variety of computer software packages that have been written specifically to solve ordinary differential equa tions both algebraically and numerically Consider then the follow ing two consecutive chemical reaction system involving three chem ical species S1 Sz and S3 v1 v2 Sy gt So gt S3 1 7 DIGITAL COMPUTER MODELS 29 We assume that each chemical reaction is governed by simple mass action kinetics so that the reaction rates v and v2 are given by v k s v k62 We can write down a simple mathematical model that describes the time evolution of the three species by invoking the conservation of mass That is dS LPNA dt 191 dS ay Ie dt 1V1 202 dS 3 ks dt 2122
28. ion currents mass flows etc In contrast to the internal variables the external variables are fixed and in principle are under the complete control of the experimenter Often the external variables are clamped to some fixed values cf voltage clamp The external variables are also called boundary variables because they are considered to be at the boundary of the system and the rest of the universe In an experiment involving liv ing cells the external variables will often be concentrations of nu trients in the medium the level of oxygen or the concentrations of externally applied drugs or growth factors Physically external variables are clamped by some kind of buffering mechanism The buffering mechanism can simply be a large exter nal reservoir so that any exchange of mass between the system and the external reservoir has a negligible effect on the external concen tration Alternatively there may be active mechanisms maintaining an external concentration A classic example of active maintenance of an external variable is the voltage clamp used in electrophysi ology Finally external concentrations may simply be slow moving compared to the timescale of the model so that over the study period the external concentrations change very little A typical example of the latter is the study of a metabolic response over a timescale that is shorter than gene expression This permits a modeler to study a metabolic pathway without considering the e
29. nts can be actually be proved to be true or false is mathematics Starting with a set of axioms mathematicians derive theorems that can be shown beyond any doubt to be true or false In contrast scientific models based on observations cannot be proved correct This is because it is simply not possible to test every pos sible circumstance in which the model may apply Instead we are left with two options 1 We falsify a model by finding one observa tion that the model fails to predict When this happens the model must be changed or abandoned For example the hypothesis RNA is never transcribed into DNA can be falsified simply by finding one instance where it happens e g the life cycle of the HIV virus Although the idea of falsifying a model is appealing in practice it is not often used Instead models are devised and used based on the degree of confidence we have in them which leads to the second option 2 We find observations that are consistent with the model In this case our confidence in the model to make useful predictions increases Finding observations that are consistent with a model or match a particular prediction is called model validation The word validate may suggest that we are suggesting that the model is cor rect but this is not true A validated model is simply one where our confidence in a model s ability to predict and provide insight has increased This is an extremely important point that is often missed The t
30. omponents to use as buffering circuits for the computational circuits described in the previous section To illustrate a simple example consider the non inverting amplifier shown in Figure 1 11 From the circuit we can determine the cur rents 7 and 72 from the voltage drops across R and R that is ur O i Vo UF l2 Ry R2 i Note that the ground voltage is 0 In an ideal op amp the amount of current flowing into positive and negative inputs is zero this means that the two currents 7 and 72 must be equal therefore it must be true that UF Vo UF UF UF VO gt 0 Ry R gt Ri R2 Another property of an op amp is that the output voltage is equal to the gain A times the difference in the input voltage vo A vr vr 1 5 ELECTRICAL MODELS 23 Rearranging this yields i vr 1 1 Inserting vr into the current equation gives ho ee e I A r a Ry R2 The final property of an op amp is that the gain A is very large ideally infinite therefore the above equation simplifies to vr vr vo am 0 Ri R2 This can now be rearranged to solve for vo 1 The remarkable thing about this result is that the gain of the op amp A has vanished Computationally this circuit will take an input volt age vz and multiply it by a scaling factor The only disadvantage to the non inverting circuit is that the scaling will be greater than one it is no pos
31. oncentrations of reactants and products show not net change over time however unlike ther modynamic equilibrium there is a net flow of mass or energy between the system and the environment At steady state the system will continually dissipate entropy to the external envi ronment while the entropy level of the system itself remains constant Transient State Under a transient state a system will be moving from either one steady state to another or from a steady state to thermodynamic equilibrium 1 11 Software for Solving ODEs There are many software tools available for solving ordinary differ ential equations We can split the list into open source and propri etary commercial tools Of the commercial tools Matlab http www mathworks com and Mathematica http www wolfram com are well knowm Matlab is a numerical computing environ ment with significant strengths in matrix manipulation and plotting For solving differential equations Matlab offers a variety of solvers including traditional Runge Kutta ode45 and ode23 and solvers specifically designed for stiff systems based on the Gear method odel5s All these solvers implement an adjustable step size al gorithm so that the execution proceeds most efficiently A more 36 CHAPTER 1 INTRODUCTION TO MODELING Thermodynamic Steady State Transient State Equilibrium Figure 1 17 Classification of System Behavior detailed discussion of different kinds of solvers and how t
32. or by analytic means Although visual models can be used to make predictions the kinds of predic tions that can be made are limited The use of mathematical models opens up whole new vistas of study which visual models simply can not match One word that has been used but has been undefined is the word System Here we define what is mean by the system System The part of the Universe we are interested in studying with a defined boundary and set of processes within Open Closed and Isolated Systems Like many other scientific studies modeling divides the world into two parts the system and the surroundings The system is that part of the universe we wish to focus our attention This might be a single enzyme reaction a pathway a single cell an organ or an entire mul ticellular organism The surroundings includes everything else The boundary between the system and the surroundings can have dif ferent properties which determines whether the system is isolated closed or open See Figure 1 2 All living systems are open be cause they exchange energy and matter with their surroundings An isolated system is one where nothing can enter or leave Such sys tem can only be approximated as it is very difficult to completely isolate a part of the Universe A thermos flask can be considered 1 2 MODELS AND SYSTEMS 5 an isolated system although an imperfect one A closed system is one that can exchange energy with the surroundin
33. pective of the mass of substance itself Thus 1 mole of glucose is the same amount as 1 mole of the enzyme glucose 6 phosphate isomerase even though the mass of each type of molecule is quite different The actual number of particles in 1 mole is defined as the number of atoms in 12 grams of carbon 12 which has been determined empirically to be 6 0221415 x 107 This definition means that 1 mole of substance will have a mass equal to the molecular weight of the substance this makes is easy to calculate the number of moles using the following relation mass moles molecular weight The concentration of a substance is expressed in moles per unit vol ume and is usually termed the molarity Thus a 1 molar solution means mol of substance in 1 litre of volume Dimensional Analysis Dimensional analysis is a simple but effective method for uncover ing mistakes when formulating kinetic models This is particularly true for concrete models where one is dealing with actual quantities and kinetic constants Conceptual models are more forgiving and don t usually require the same level of attention because they tend to be simpler Amounts of substance is usually expressed in moles and concen trations in moles per unit volume mol Reaction rates can be expressed either in concentrations or amounts per unit time depend ing on the context mol t mol t7 Rate constants are expressed in differing units depending on the form of
34. r for calculations predictions etc This definition embodies a number of critical features that defines a model the most important is that a model represents an idealized description a simplification of a real world process This may at first appear to be a weakness but simplification is usually done on purpose Simplification is important because it allows us to com prehend the essential features of a complex process without being burdened and overwhelmed by unnecessary detail In the water tank model we assumed for example that as the water flowed from each tank any resulting water turbulence had only a negligible effect on the system dynamics Models come in various forms including verbal written text visual mathematical and others Biology particular molecular biology has 4 CHAPTER 1 INTRODUCTION TO MODELING a long tradition of using visual models to represent cellular struc ture and function one need only look through a modern textbook to see instances of visual models on every page Visual models have been immensely useful at describing complicated biological processes but are limited in their scope A more interesting way to describe models is to use mathematics a language designed for logical reasoning Mathematical models are useful in biology for a number of reasons but the three most important are increased pre cision prediction and the capacity for analysis Analysis can be carried out either by simulation
35. rinsitic stochastic there are additional properties of models that can be used to categorize them Table 1 3 1 Linear or Non Linear 2 Dynamic or Static 3 Time invariant or time dependent 4 Lumped or distributed parameter models Table 1 3 Additional categories of models Dynamic and Static Models A static model is one where the variables of the system do not change in time For example a circuit made up of only resistors can be modeled as a static system because there are no elements in the circuit that can store or dissipate charge thus currents and volt ages are considered instantaneous without any time evolution Most interesting model are dynamic 1 4 CLASSIFICATION OF MODELS 15 Time Invariant Systems A time invariant model is one where the model does not explicitly depend on time This means that given a time invariant model run ning the model at 0 or tf 10 makes no different to the time evolution of the model If a parameter of the system depends on time then the model is called time dependent An example of a time dependent model is where we apply a drug in the form of a pulse and the duration of the pulse depends on when the drug is admin istered An example of a time dependent non biological model is a parking lot where the price of a ticket depends on the time of day We will have more to say about time invariant systems in a later chapter when we will talk about linear time invariant systems LTD
36. s either negligible or unimportant In many non biological systems such an approximation might be quite reasonable However biological sys tems are particulate in nature and operate at the molecular level As a result biological systems are susceptible to noise generated from thermal effects as a result of molecular collisions For many systems the large number of particles ensures that the noise generated in this way is insignificant and in many cases can be safely ignored For some systems such as prokaryotic organisms the number of par ticles can be very small In such cases the effect of noise can be significant and therefore must be included as part of the model 1 10 Behavior We can classify the behavior of systems in to three broad categories thermodynamic equilibrium steady state or transient The steady state can be further classified into stable or unstable steady states For example if we pick biochemical systems as our example then 1 11 SOFTWARE FOR SOLVING ODES 35 Thermodynamic Equilibrium In this state the concentrations of reactants and products show not net change over time in ad dition the rates of all forward and reverse reactions are equal This means that at thermodynamic equilibrium there is no net movement of mass from one part of the system to another and not net dissipation of energy In thermodynamics equilibrium is the state that maximizes the entropy of the system Steady State At steady state the c
37. sible to scale vz by a factor lt 1 More importantly given the high input and low output impedance the circuit doesn t suffer the loading issues of a simple resistive cir cuit Combining op amps with resitive capacitive circuits is a com mon way to make analog computational modules that can be easily connected together without the modules interfering with each other Equation 1 1 is of interest because it helps to understand the opera tion of the circuit intuitively Given that the gain A of the op amp is huge we can approximate equation 1 1 to Ul UF 24 CHAPTER 1 INTRODUCTION TO MODELING Figure 1 11 Non Inverting Op Amp Circuit This tells us that what ever the input voltage vz is vr will match it In order for this to happen as the input voltage changes vz the op amp must be going through some some kind of dynamic change so that vp will match any change in vz The explanation is simple but crucial to understanding the operation of the circuit First we must recall that an op amp is a difference amplifier that is Vo A vy v_ or vo A v vp Imagine we increase vz above vpr this means the difference vy vf increases which in turn is amplified resulting in an increase in vg If vo increases then part of this is passed back via the resister R2 to vp so that vr increases If vp increases the difference vy vr decreases which in turn reduces the rise in vg The feedback cycle continues until
38. ssed as a ratio of lengths which will by necessity have the same dimension 1 9 Approximations By their very nature models involve making assumptions and ap proximations The best modelers are those who can make the most shrewd and reasonable approximations without compromising a mod el s usefulness There are however some kinds of approximations which are useful in most problems these include 1 9 APPROXIMATIONS 33 e Neglecting small effects e Assuming that the system environment is unchanged by the system itself e Replacing complex subsystems with lumped or aggregate laws e Assuming simple linear cause effect relationships where pos sible e Assuming that the physical characteristics of the system do not change with time e Neglecting noise and uncertainty Neglecting small effects This is the most common approximation to make In many studies there will always be parts of the system that have a negligible effect on the properties of the system at least during the period of study For example the rotation of the earth the cycle of the moon or the rising and setting of the sun will most likely have negligible influence on our system Assuming of course we are not studying circadian rhythms Assuming that the system environment is unchanged by the sys tem itself This is a basic assumption in any study The minute a system starts to affect the environment we have effectively extended the system boundaries to
39. te out the differential equations for this system c Assign arbitrary but realistic values to the various param eters in the model and compute the steady state that is when d dt 0 Report the steady state concentrations of S4 and S2 d What is the effect of doubling all E by a factor on the steady state concentrations of S and S2 Pharmokinetic Models Figure 1 22 shows a simple compartmental model of what happens to a drug when it is injected intravenously into a pa tient The model comprises of three compartments the blood stream the liver and the actual site of action Itis assumed that the drug can freely exchange between the liver and blood compartments but that the liver can also irreversible degrade the drug In addition some drug is lost by excretion through the kidneys Finally some drug accumulates irreversibly at the site of action 44 CHAPTER 1 INTRODUCTION TO MODELING i Concentration Time Vy Input Profile at this Point Accumulation Exchange Vs 3 v4 lt _ _ _ _ _ lt rr v2 Degradation Us Excretion Figure 1 22 Simple model of a drug begin injected and dis tributed through the body We will assume that movement of the drug between compart ments obeys simple first order kinetics For example the rate v2 at which drug enters the liver is given by k2 Dp where Dy is the concentration of drug in the blood stream and k3 is the kinetic constant for the process
40. tta method See Chapter gt gt x0 8 8 27 gt gt tspan 0 20 gt gt t x ode45 lorenz tspan x0 gt gt plot t x The first line sets up the initial conditions for the three variables the second up the duration of the simulation and the third line carries out the actual simulation The last line plots the data as a graph The next example shows how we can use Python to solve the Lorenz system This relies on the SciPy and NumPy packages The SciPy package includes a wide variety of numerical methods for solving different kinds of problem The first part of the script imports the required libraries import scipy from scipy integrate import odeint import numpy import matplotlib pyplot as plt The second part of the script defines the system of differential equa tions def lorenz x t sigma 10 rho 28 beta 8 0 3 return sigma x 1 x 0 x 0 rho x 2 x 1 x 0 x 1 beta x 2 Note the indentation of the lines this is important for Python Fi nally we invoke the ODE solver and plot the resulting simulation data 40 CHAPTER 1 INTRODUCTION TO MODELING xInitial 0 1 1 05 t numpy arange 0 30 0 01 lorenzSolution scipy integrate odeint lorenz xInitial t plt plot lorenzSolution plt show Domain Specific Tools In addition to tools that can solve differential equations there are also domain specific tools that allow a user to specify the phys
41. vp converges to vz at which point the circuit reaches a steady state that is vo stops changing and vy vF Assignment 2 Using the same approach to what was used in the non inverting amplifier Figure 1 12 derive the relationship between the input and output voltage for the inverting amplifier 1 Show that the inverting amplifier behaves like a voltage divider while at the same time inverting the signal 1 5 ELECTRICAL MODELS 25 2 Show that the voltage at vp is zero Using the explanation that vy uF for the non inverting op amp explain why vp is zero in an inverting op amp Figure 1 12 Non Inverting Op Amp Circuit Figure 1 13 shows three computational op amp circuits a summer an integrator and a differentiator Vo Summer Integrator Differentiator Figure 1 13 Three different computational circuits made from op amps Note the similarity with the non op amp circuits in Figure 1 7 26 CHAPTER 1 INTRODUCTION TO MODELING 1 6 Block Diagrams In may text books block diagrams are used to illustrate a model We already seen that electrical circuits can be used perform vari ous computations and these can be represented in the form of block symbols Figure 1 14 shows some common block symbols found in the literature Using these symbols let us draw a block diagram that represents the differential equation and its solution Solving a differential equation involves integration Therefore the first
Download Pdf Manuals
Related Search