Equation Help, Examples & Guidelines

-1-

TABLE OF CONTENTS 1

EQUATION HELP ......................................................................................................................................... 3 1.1 Equation Types ....................................................................................................................................... 3 1.2 Built-in Functions.................................................................................................................................... 5 1.3 Operators................................................................................................................................................. 7 1.4 Hints........................................................................................................................................................ 7 2 EXAMPLES.................................................................................................................................................... 8 2.1 5-Minute Tutorial .................................................................................................................................... 8 2.2 Calcium Oscillations - Phase Plane Plots.............................................................................................. 10 2.3 Cardiac (Purkinje fiber) Action Potentials - Stiff Equations................................................................. 13 2.4 Cell Permeability - Curve fitting Empirical Data.................................................................................. 16 2.5 Chemical Kinetics Module.................................................................................................................... 19 2.6 Hodgkin Huxley Model of Nerve Axon - Submodels........................................................................... 22 2.7 Oxygen Diffusion from Blood Capillary - Boundary Value Problem (Shooting Method) ................... 24 2.8 Local Anesthetic - Oscilloscope Mode ................................................................................................. 26 2.9 Using indexed variables in a linear array of cells.................................................................................. 28 2.10 Monte Carlo Model of a single Na channel channel under a patch clamp ............................................ 30 3 MATHEMATICAL EXAMPLES ................................................................................................................ 32 3.1 Custom DT Method............................................................................................................................... 32 3.2 Custom DT Method............................................................................................................................... 33 3.3 Difference Equations............................................................................................................................. 34 3.4 Floating Point Exceptions ..................................................................................................................... 35 3.5 NETFLOW Built-in .............................................................................................................................. 37 3.6 Phase Plane ........................................................................................................................................... 38 3.7 Relations and Logical Operators ........................................................................................................... 39 3.8 Stiff Systems ......................................................................................................................................... 40 3.9 Vector Finals ......................................................................................................................................... 41 4 HOW DO I .................................................................................................................................................... 42 4.1 Control Plotted Variables ...................................................................................................................... 42 4.2 Define and Use Sliders.......................................................................................................................... 44 4.3 Display Numerical Values .................................................................................................................... 45 4.4 Make a Phase Plane Plot ....................................................................................................................... 46 4.5 Perform a Curve Fit............................................................................................................................... 47 4.6 Plot a Family of Curves......................................................................................................................... 51 4.7 Solve Boundary Value Problems .......................................................................................................... 52 4.8 Chemical Kinetics Module.................................................................................................................... 55 4.9 Using Conveyors................................................................................................................................... 58 4.10 Use Indexed Variables (Arrays) ............................................................................................................ 59 4.11 Using Multidimensional Arrays ............................................................................................................ 60 4.12 Using Ovens.......................................................................................................................................... 61 4.13 Using Queues ........................................................................................................................................ 63 4.14 Using the Flow Chart ............................................................................................................................ 65

-2-

1 1.1

EQUATION HELP EQUATION TYPES

Integration Method METHOD name

name = Euler, RK2, RK4, Auto, or Stiff

Integration Parameters (Default Settings) STARTTIME = 0 STOPTIME = 10 DT = 0.02 DTMIN = 1.0e-6 DTMAX = 1.0 TOLERANCE = 0.01 DTOUT = 0 ROOTTOL = 0.001

Initial time Final time Integration time step (Euler, RK2, and RK4 methods only) Minimum DT (Auto and Stiff methods only) Maximum DT (Auto and Stiff methods only) Relative accuracy for Auto and Stiff integration methods Output time interval (0 = store every step) Relative accuracy for built-in root finders

You can change the names of the above symbols and TIME by using the RENAME statement like this: RENAME TIME = X RENAME STARTTIME = Xo RENAME STOPTIME = Xf

Initialization Equations x(STARTTIME) = expression INIT x = expression INIT(x) = expression

Differential Equations x' = expression d/dt(x) = expression FLOW x = expression x(t) = x(t - dt) + (expression) * dt x = x + dt * (expression) You can define second-order (and higher) equations directly using the "prime" notation like this: x'' = expression Berkeley Madonna internally translates this into a system of two first-order equations with variables named x' and x. You can refer to these variables in equations just like other variables you define. You must provide initializers for each of these variables. For example, here is a third-order system: x''' init init init

= -2*x'' - 3*x' - x + 1 x'' = 2 x' = -1 x = 0

Difference Equations x(t + dt) = expression NEXT x = expression The value of x is replaced by the value of the expression on the right-hand side. Therefore, to implement a typical finite difference equation where some value f(x) is added to the current value of x, you must write: x(t+dt) = x + f(x) NEXT x = x + f(x)

or

-3-

Limit Equations LIMIT v >= c LIMIT v = c Constrains variable v to be = ... = t THEN h ELSE 0 PULSE(v, f, r) Impulses of volume v, first time f, repeat interval r SQUAREPULSE(t, d) Pulse of height 1 starting at time t with duration d GRAPH(x) (x1,y1) (x2,y2) (x3,y3) ... Lookup x in list of points using linear interpolation

Array Functions ARRAYSUM(x[*]) ARRAYMEAN(x[*]) ARRAYSTDDEV(x[*])

Sum of all elements in array x Average of all elements in array x Standard deviation of all elements in array x

-6-

Dataset Functions Any dataset in your model can be used as a piecewise-linear function. Simply use the dataset as if it were a built-in function taking one or two arguments: t = #temperature(x) q = #charge(x,y) This example assumes #temperature is a vector (one-dimensional) dataset and #charge is a matrix (two-dimensional) dataset. Berkeley Madonna computes the result by performing linear interpolation on the points nearest the argument value(s).

1.3

OPERATORS

Arithmetic Operators x + y x – y x * y x / y x ^ y x ** y -x

Addition Subtraction Multiplication Division Exponentiation Same as x ^ y Negation

Relational Operators x x x x x x

= y y < y y >= y

Equal to Not equal to Less than Less than or equal to Greater than Greater than or equal to

Logical Operators x AND y x OR y NOT x

Logical AND Logical OR Logical inverse

Conditional Operator IF x THEN y ELSE z

1.4

HINTS

Making Non-negative Reservoirs and Unidirectional Flows Use LIMIT statements to make reservoirs non-negative and flows unidirectional. example: r(t) = r(t - dt) + (i - o) * dt LIMIT r >= 0 LIMIT i >= 0 LIMIT o >= 0

Makes reservoir r non-negative Makes inflow i uni-directional Makes outflow o uni-directional

Creating Periodic Functions y = F(MOD(TIME, i))

F repeats over the time interval i

Example: Square wave with period = 1 y = IF MOD(TIME, 1) >= 0.5 THEN 1 ELSE 0

-7-

For

2 2.1

EXAMPLES 5-Minute Tutorial

Simple Chemical Reaction: A → B → C 1. Enter (Type or drag & drop or copy & paste) the following equations into the Equation Window (under the Model menu, or Ctrl-E) METHOD RK4 STARTTIME = 0 STOPTIME = 10 DT = 0.02 d/dt (A) = - ka*A d/dt (B) =ka*A - kb*B d/dt (C) = kb*B init A =100 init B = 0 init C = 0 ka =1 kb =1 2. Click the Run button, or type Ctrl-R Madonna automatically plots the first two variables in the Equation window. The first is scaled on the left, the second on the right. The next six variables are easily addressed by buttons on the left bottom margin. 3. Use the buttons on the left bottom margin of the display to toggle the graphs of A,B, and C. Toggle C on. Use shift-click B to change the scaing axis for B from right to left axis. General control of ploting variables (e.g. adding other variables) is handled by a dialog initiated after double clicking anywhere in the x-y plane. 4. The Graphic tool bar is located at the top of the graph window.

Click the button to see the graph change into a table of numerical values. Toggle the legend. See appendix below (step 9) for a general description of all buttons.

button for a

5. Open the Parameter Window (Parameters Menu) using the menu or Ctrl-Shift-P:

Change the value of ka by double-clicking it and entering the changed value into the edit box next to the Reset button. Parameters that have been changed are indicated with an asterisk - they can be reset to default values by the Reset button. 6. Open the Define Sliders dialog (in the Parameters menu). Double click on ka and on kb to add them to the Slider list. Click OK. Move the sliders back and forth: when you release the mouse, Madonna automatically runs the program with the current parameter values. Double click within the text area of the sliders to open the Sliders dialog where you can add sliders, adjust the current slider limits etc. -8-

7. Plot a family of curves. Click on Batch Runs in the Parameters menu and set it up as shown. In this case we want to make several automatic runs for different values of ka, so we select ka in the Parameter popup. We will display 10 runs with values ranging from ka = 0.1 to 1.0. The intermediate values are automatically calculated and displayed in the Values: scroll window.

Click OK. The concentration profiles are calculated and displayed. Toggle the L button and the ka corresponding to each profile is displayed in the legend. Toggle A and C off the graph so that you can see the unencumbered behavior of B. The same can be done for A and C.

8. Parameter plots. Automatically plot concentration "behavior" as a function of the parameters. Open the Parameter Plot dialog (under the Parameters menu). Select ka from the Parameter: list (at the top of the dialog). Set # of Runs = 20, initial value = 0.1, final value = 5.0. From the Variables list, choose B by double clicking on it. Check the Maximum. box (lower right). Click on Run. Toggle the tool bar button with the dot in the middle to show the actual points that were used to calculate the curve. 9. APPENDIX

-9-

2.2

Calcium Oscillations - Phase Plane Plots

Features: • Limit Cycles • Phase Plane Plots Description of the Model: When stimulated by hormones or neurotransmitters, many types of cells burst into repetitive spikes of intracellular calcium release. The period of these oscillations ranges from less than 1 second to more than 30 minutes. These oscillations are thought to be an important attribute of intra and intercellular signaling mechanisms. The following model was taken from: Goldbetter, A., G. DuPont, M. Berridge (1990). Minimal model for signal-induced Ca2+ oscillations and for their frequency encoding through protein phosphorylation. Proc. Natl. Acad. Sci. USA 87:1461-65. In the calcium signaling system, shown below, a ligand binds to a receptor, R, which, after a series of steps, produces intracellular inositol triphosphate (IP3). IP3 then triggers the release of calcium (Ca++) from an intracellular storage vesicle. However, not all Ca++ containing vesicles are controlled by IP3. This model considers the case where the concentration of IP3 is constant (steady state) and shows that the interplay of Ca++ sources and sinks can lead to limit cycle oscillations.

Letting fl denote the normalized concentration of IP3, we write conservation equations for the concentration of intracellular calcium, Z, and the concentration in the IP3-insensitive pool (pool 2), Y:

The fluxes into and out of the IP3 insensitive pool (2) are the key nonlinearities controlling the behavior of the system.

- 10 -

Follow These Steps: 1. Run the program 2. Replot the results in a phase plane a. Double click anywhere in the x-y plane of the graph (or select Choose Variables from the Graph menu.) The following dialog appears...

b. Remove Z from the Y axis. c. Replace TIME on the x axis with Z by selecting Z in the drop down X Axis menu. d. Click OK. You will see a phase plane plot starting from the initial condition and moving into a closed loop limit cycle. 3. Change initial conditions. a. Toggle the (overlay) button on the graph tool bar. (This allows results to be compared, but it is not required.) b. Toggle the (initial conditions) button. The graph cursor changes from to . c. Place the cursor anywhere in the Y-Z plane and click the mouse. Madonna interprets this as a new set of initial conditions for Y and Z given by the coordinates of the point that you have selected. Madonna automatically runs the program with these new conditions; numerical values are displayed in the graph window, just above the graph. d. Repeat this (c.) at several different positions. In each case you will see the plot converge on the same limit cycle. 4. Explore a parameter range for oscillations. a. Toggle the button to turn overlay off. b. Select Batch Runs ... from the Parameters menu. Fill in the dialog (below) as shown: Set Parameter = beta, # of Runs =10, Initial Value = 0.1, Final Value = 1.0 c. Click OK. You will find that only values of beta (IP3) ranging from 0.3 to 0.7 give stable limit cycles. (some of the other plots may not move much or they may be buried)

- 11 -

- 12 -

2.3

Cardiac (Purkinje fiber) Action Potentials - Stiff Equations

Feature: • Rosenbrock method for stiff equations • Oscilloscope Mode • Saving memory with large detailed calculations Description of the model: A normal heart beating at say 75 beats per minute generates spontaneous action potentials every 0.8 sec. During the major portion of this period the membrane potential, Vm, changes slowly, but the beginning of each action potential is characterized by a precipitous large rise in Vm which is completed within a msec. As a result cardiac action potential models are generally stiff and ill conditioned. We illustrate use of the Stiff solver with the DiFrancesco - Noble (1985) model. It consists of 16 highly non linear differential equations with coefficients that are complex transcendental functions of Vm. In addition we will consider a case where an anti-arrhythmic drug (e.g. lidocaine) is applied to the heart beginning at time =0. (Starmer et. al. 1984). The model simulates how the action potential progressively changes as the drug, binds to and blocks the sodium channel. Note that the time scale in this problem is in seconds (contrasted with nerve excitation where the time scale is typically msec) Follow these steps: 1. Run the model. Use the zoom to see more detail on any portion of the plot. Return to the default scaling by toggling the z button on the display tool bar. (The Z toggle will take you back one zoom; if you zoom twice in a row you will have to toggle Z twice in a row to return) 2. Superimpose the results. Use the T (trigger) button on the tool bar to toggle in and out of oscilloscope mode which is set in this model to trigger at Vm > -60) 3. Check the integration parameters. Activate the Parameter Window by selecting it in the Parameters menu. The default integration algorithm used in this simulation, Rosenbrock (stiff), is indicated at the top of the window. This algorithm continually adjusts the step size (dt) to optimize the computation speed within a specified degree of accuracy. When variables change fast, dt becomes small, when they change slow, dt is large. The method requires 3 computation parameters that are not relevant with more common fixed step size (Euler, Runge-Kutta) methods. These are: a. DTMIN = the smallest dt that the algorithm will use b. DTMAX = the largest dt that the algorithm will use c. TOLERANCE = the specified error limit (accuracy) (see the manual for further details)

4. Change the integration algorithm. To test the effectiveness of the Rosenbrock method for this computation, we repeat it with the more conventional Runge Kutta. Bring up the Parameter Window , and use the popup menu at the top to change from Rosenbrock (stiff) to Runge-Kutta4

- 13 -

5. Use DTOUT to prevent memory problems. If you try to run the program as is you may encounter an out of memory error message. The problem is that in default mode (DTOUT = 0) Madonna tries to store every computational point and we are set to store STOPTIME/DT = over a million points. You can control the storage by using DTOUT to specify the interval between stored points. If, for example, you set DTOUT = 0.1, Madonna will continue to compute with dt = 5e-5, but it will only store 800 points at 0.1 sec apart for plotting. Now run the program and note how long it takes! (Any substantial increase in dt generates an overflow). In addition, if you used a non aero DTOUT, you may find that the plot has been severely compromised. Zoom in on a peak to verify this. 6. Plot dt. To inspect how the stiff algorithm changes each dt, change the integration method back to Rosenbrock (stiff) and change STOPTIME (Parameters Window) to 2 so that details of a single action potential will be apparent. Run the model and, toggle the dt_ button on the display and values of dt (scaled on the right axis) will be superimposed on the action potential. The variable dt_ was defined within the Equation Window by the statement dt_ = stepsize where stepsize is a reserved term that makes the program dt available for user manipulation.

Model Equations Consult the original articles for details: The cardiac model was developed for rabbit purkinje fibers by DiFrancesco, D. and D. Noble (1985). "A model of cardiac electrical activity incorporating ionic pumps and concentration changes." Philos Trans R Soc Lond B Biol Sci 307(1133): 353-98. Corrected Equations were kindly provided by C. Cabo and R. C. Barr (personal communication) - also see: Cabo, C. and R. C. Barr (1992). "Propagation model using the DiFrancesco-Noble equations. Comparison to reported experimental results." Med Biol Eng Comput 30(3): 292-302. The anesthetic model was developed for Na channels in squid axons by Starmer, CF; Grant, AO; Strauss, HC. "Mechanisms of use-dependent block of sodium channels in excitable membranes by local anesthetic" Biophysical Journal, 1984 Jul, 46(1):15-27.

-------------------------------------------------------------------------The tutorials are designed to introduce many of Berkeley Madonna's features by leading you step by step through a computation. In many instances the examples will be taken from our laboratory course, which can be accessed on the web at: http://mcb.berkeley.edu/courses/mcb137/.

- 14 -

Modeling the Cell Cycle The protein cyclin is a key ingredient in the cell cycle. Its periodic buildup and breakdown in the cell parallels the cycle, and may actually drive it. When cyclin exceeds a certain threshold it begins to combine with and activate a protein kinase (called CDC2 kinase, not shown) to form a complex called "maturation promoting factor" (MPF), which appears to stimulate mitosis. The CDC2 kinase stimulates degradation of cyclin by activating a protease. Thus the system has a negative feedback: cyclin stimulates its own degradation via MPF and protease. The model presented here is based on: Goldbeter, A. (1991). A minimal cascade model for the mitotic oscillator involving cyclin and cdc2 kinase. Proc Natl Acad Sci. 88:9107-11. Denote the concentration of cyclin by C, MPF by M, and protease by X. A diagram of the kinetic scheme is shown below.

The model equations are:

- 15 -

2.4

Cell Permeability - Curve fitting Empirical Data

Features: • Importing empirical (numerical) data • Curve fitting model parameters to the data. Description of the model: Experimental data was obtained by measuring how cell volume, v, responds to the sudden introduction of a permeable solute into the bathing medium. The problem is to extract both the water permeability (Pw) and the solute permeability (Ps) from the data. The experiment is illlustrated in the figure. The model assumes that all movements are rate limited by moving across the membrane barrier; it takes account of osmotic water flow coupled to solute diffusion through the membrane. Equations are described below.

Follow these steps: 1. Run the model. Curve fitting requires an initial guess of all unknown parameters. In our case we guess Pw = 0.002 and Ps =0.001. These are entered in the model as though they were known (see Equation Window). Running the model shows the results of our guess verifies that the solution shows first a volume, v, shrinkage and then a swelling back towards normal. 2. Import data a. Select Import Dataset from the File menu. The standard oOpen dialog will appear b. Open the enclosed data file entitled "v_Data". This is a simple 2 column text file (Data can be text or spread sheet files with an unlimited number of columns). The dialog shown below will appear. It describes the data file, allows you to rename it, and allows you to choose which column represents the independent (x axis) variable and which represents y. Numerical data files in Berkeley Madonna are referred by by a # prefix. The default name is shown below. c. Click OK

- 16 -

Notice the large discrepancy between the data and the model run with the guessed parameters (note change in scale). 3. Curve fit a. Click Curve fit in the Parameters menu. The dialog shown below (truncated on the right side) appears. b. Double click on the unknown parameters Pw and Ps. c. Use the FIt Variable pulldown menu to select v. d. Check that #v_Data is entered in To Dataset. e. Click OK. You are finished! Parameter values used in the plotted curve fit can be found in the current Parameter window. (See the Users manual for more options on curve fitting. For example, the right hand side of the dialog can be used to limit the permissible parameter values.)

Model Equations The model is described by two coupled differential equations: d/dt(nw) = Pw*A*(Osmi-Osmo) d/dt(ns) = Ps*A*(Co - Ci) where: nw = number of moles of water inside the cell Pw = molar water permeability of cell A = area of cell Osmi = osmolarity inside cell Osmo = osmolarity outside cell ns = number of moles of permeable solute inside the cell Ps = solute permeability Co = permable solute concentration outside Ci = permable solute concentration inside

- 17 -

The data was obtained from red blood cells made into "ghosts" where the internal solutes were replaced by saline making the cell volume virtually equal to the internal water volume. It follows that Vol = nw* Vw Vol = cell volume in cm^3 Vw = molar volume of water v = Vol*1.0e12 = cell volume in µ^3 The rest of the model consists of initial conditions, definitions, and conservation relations. The size of the model could be reduced through algebraic substitutions; but why bother? Let Madonna figure it out!

- 18 -

2.5

Chemical Kinetics Module

Transient State Michaelis Menten Kinetics Features: The chemical reactions module is a template that can be appended to Berkeley Madonna models. It allows you to enter reactions in ordinary chemical notation; Berkeley Madonna will automatically set up the corresponding kinetic equations base on the law of mass action and enter them into the equation window. The equation window is fully accesible for writing additional equations. Description of the model: Non steady-state Michaelis Menten model for enzyme kinetics is set up and integrated. Follow these steps: 1. Open the Chemical Reactions dialog by selecting it from the Model menu: Modules popup. 2. Using ordinary chemical notation, enter the first reaction in the Reactants and Products boxes as shown below. 3. Enter Kf, the forward rate constant, and Kr, the reverse rate constant, as shown. 4. Click Add, and the reaction will be recorded in the Reactions list. 5. Enter the Initial Concentrations as shown.

6. Repeat the same operation for the second reaction:

- 19 -

7. Click OK. The reactions are displayed in the equation window:

The reaction objects (gray areas) display Berkeley Madonna's interpretation of the Chemical Reactions dialog. Change the STOPTIME to 300. (All parameters can be changed in the usual way with the Parameters window.) Click the Run button to see the results. It is instructive to plot the results against log (time). Double click on the x axis to open a dialog for making this change. Important Notes: Contents of the reaction objects can be supressed by the Show check boxes on the lower right in the Chemical Reactions dialog. Try them. Edit the reactions. Double click on a reaction object to reopen the reactions dialog. The reaction you clicked is selected. Change it in the Reactions and Products boxes at the top. After a reaction is edited, you must click the Modify button for the change to take effect. If the Add button is used Berkeley Madonna will interpret it as an additional reaction. To remove a selected reaction use the Remove button.

- 20 -

To change the mass action interpretation of any reaction rate simply add your own rate equation in the Equation window. E.g. if the desired first reaction rate is given by f(S,E), then write RXN1 = f(S,E) in the Equation window. As long as your equation comes after the reaction object, it will supercede the equation in the object. (If your equation comes before the reaction object, your definition will be superceded by the one supplied in the object.)

- 21 -

2.6

Hodgkin Huxley Model of Nerve Axon - Submodels

Features: • Flowchart interface with serious model • Use of sub-models Description of the model: This is the classical model of excitation that established the primary paradigm for cellular neuurobiology for the past 60 odd years. It describes the excitation properties of nerve axons in terms of ionic movements across its membrane. (Hodgkin, A. L. and Huxley, A. F (1952). "A quantitative description of membrane current and its application to conduction and excitation in nerve." J. Physiol., 117: 500-544.) Follow these steps: 1. Run the model to see the development of an action potential. 2. Change the stimulus intensity by moving the slider. When it is below a threshold value, the action potential disappears leaving only the stimulus "artifact". Follow the model Flowchart 1. The top flow chart shows the flow of charge, Q, to the inner membrane surface. It is composed of a balance of stimulating current counted positive in the inward direction -i.e. into the axon) and ionic currents carried by Na, K, and "leakage" ions (counted positive in the outward direction) The charge is related to voltage (membrane potential E) through the membrane capacitance. 2. Sub model icons INa_, IK_, and IL_, are used to represent the ionic currents. Double click on the IK_ icon to locate IK. It is set proportional to the difference between two driving forces, the membrane potential E, and the K equilibrium potential, EK. (The equilibrium potential reflects the free energy contribution of the K concentration gradient.) IK = gK (E-EK) where gK is the K conductance of the membrane. gK is not constant; K leaves the cell through channels whose status (open or closed )is controlled by 4 identical gates (called n gates) that are voltage dependent.. All 4 gates must be open for the channel to conduct. Letting n denote the probability of an open state for one gate, then the probability of finding all 4 gates open is n^4, and gK is given by gK = gKmax *n^4 where gKmax is the K conductance when all channels are open. 3. The n gates sub model is used to compute n. Double click on the n gates icon. You will see that n is calculated by simple first order kinetics with the caveat that the "production" and "decay" rate constants, an and bn, depend on E. The explicit dependence used in the model was determined empirically from voltage clamp data. 4. The formal treatment of INa is very similar with the exception that one of it's gates, the h gate, is basically different than the other three. This can be seen by comparing the expressions for the rate constants ah and bh with the corresponding rate constants, am and bm, for the m gates. Model Equations For a detailed exposition of the equations see any of the following: Hille, B. (1992) Ionic channels of excitable membranes. 2nd ed. Sunderland, Mass. : Sinauer Associates, 1992. Junge, D. (1992) Nerve and muscle excitation. 3rd ed. Sunderland, Mass. : Sinauer Associates, c1992. Hodgkin, A. L. and Huxley, A. F (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 117: 500-544.)

- 22 -

- 23 -

2.7

Oxygen Diffusion from Blood Capillary - Boundary Value Problem (Shooting Method)

Features: • Boundary value problem • Insertion of higher order differential equations directly into Berkeley Madonna Description of the model:

Oxygen diffuses radially from the surface of a capillary where its concentration is held constant at 100 units into tissue that consumes it at a constant rate Q units/sec. The capillary has a radius r = 4 µm. Diffusion within the the capillary is ignored, as is the variation of oxygen concentration along the z axis (These conditions follow if the flow of blood through the capillary is rapid) The flux is zero at the boundary of the tissue: r = R = 12µm,. This is the classical Krogh capillary model - it is oversimplified (e.g.blood flow is not that rapid). Run the Model Note that the boundary condition dC/dr = C' = 0 at r = 12 is satisfied. (C' is scaled on the right axis.) Model Equations The steady state diffusion equation in cylindrical coordinates, with sink Q, is given by

Higher order derivatives can be processed by Madonna if they are written with "prime" notation. (Also , within the equation window we change the name of the independent variable from the default TIME to r by the statement: RENAME TIME = r ). In "prime" notation the diffusion equation becomes

Enter this into Madonna's equation window in the form:

The integration is carried out from the capillary surface at r = 4 to the tissue boundary at r = 12. Boundary Conditions are:

The first boundary condition is also an initial condition for C. (Note that the initial boundary corresponds to r = 4 which is entered into the equation window by setting STARTTIME =4) Setting initial and boundary values: The initial condition for C' is not known. It is arbitrarily set at -2. The BVP solver will look for the appropriate initial condition; it will keep changing this condition and running the model until it finds one where the flux = 0 at r = 12. This is activated by the Boundary Value ODE dialog (Model menu Modules popup). The setup in the boundary value dialog box is shown below. To set the boundary condition: - 24 -

1. Select C' from the Set menu in the drop-down menu at the top of the dialog 2. Fill in the 2 edit boxes that lie immediately below ("=" and "at x") 3. Click on the Add button. This will automatically enter the condition in the Boundary Conditions list. 4. At the bottom of the dialog, choose which parameter to vary in order to satisfy the boundary condition. In our case, we select INIT C' since we have assigned an arbitrary value to it. 5.Click on the Add button. This will automatically enter the parameter in the Unknowns list. 6. The Minimum and Maximum edit boxes can be adjusted to narrow the search by setting upper and lower limits to the unknown.

- 25 -

2.8

Local Anesthetic - Oscilloscope Mode

Features: 1. Demonstrates the display operating in both the "Y-T recorder" and "Oscilloscope" modes. 2. Uses the mod function to generate periodic square pulses Model A local anesthetic (e.g. lidocaine) exerts its action by blocking Na+ channels in excitable cell membranes. The anesthetic can only access the blockage site when the channel is in certain "open" configurations associated with activity. This confers a useful specificity to the anesthetic; it selectively blocks the most rapidly firing cells (e.g. sensory axons mediating pain, and runaway arrhythmic cells in the heart). The model simulates an experiment where the anesthetic (acting over a time scale of seconds) is applied to a nerve axon. Progress of the anesthetic is assayed every 0.1 second with a voltage clamp. The clamp measures current passing through the Na+ channel following a jump in membrane potential from Erest = -100 mv tp Eset = 0 mv. Meaningful measurements are complete within the first 5 msec of each 0.1 sec interval because the normal channel cannot be sustained in an open configuration for more than a few msec. Use of the oscilloscope plotting mode allows the series of 5 msec measurements to be superimposed despite the fact that the entire simulation lasts 1000 msec Follow these steps 1. Run the model. The channel currents show up as negative spikes every 100 msec. 2. Toggle the (trigger) button on the top of the graph. The display now changes; the first 2 msec of each 100 msec interval are now superimposed. How Oscilloscope Mode is setup 1. Set the trigger Oscilloscope Mode causes the graph to start drawing at zero on the X-axis whenever a trigger event occurs. We want a trigger event each time we test the Na channel by jumping E from Erest (= -100) to Eset (= 0). To implement this we simply place the following statement anyplace within the Equation Window: Trigger = E > Erest and a "trigger event" will occur each time E > Erest. 2. Set the duration of the oscilloscope plot. The "Y-T recorder is set to plot from time = 0 to time = STOPTIME (1000 msec). The "Oscilloscope" has its own independent time sweep and this includes the duration of each sweep. This duration was set to 2 msec by double clicking on the time axis while the oscilloscope display was active, and then filling in the desired limits for the x- axis in the ensuing dialog. (This action requires that the Auto scale check box (arrow in the figure) is toggled to the unchecked, off position))

Generating periodic Functions To simulate the voltage clamp membrane potential that rises to 0 for the first 5 msec of each 100 msec interval and then falls to -100 mv, we first setup a sawtooth with period 100 and slope 1. Sawtooth = MOD(Time,1) Use the T button to return to the Y-T recorder and then toggle the Sawtooth variable on to verify that it behaves as expected. (Use the zoom to see more detail)

- 26 -

Next we use the sawtooth as the periodic "time base " for our clamp, i.e. E = IF Sawtooth < 5 THEN 0 ELSE -100 Note: the sawtooth, the trigger, and the duration of the oscilloscpe plot are all independent of each other. Model Equations The Hodgkin-Huxley Equations simulate a voltage clamp experiment. Action of the anesthetic is assumed to be rate limited by binding to the site which is modeled with conventional mass action assumptions taking account of the availability of the open channel configurations at all times. Blockage is assumed to be all-or-none. For details see: Starmer C. F. Grant A. O. and Strauss H. C. (1984), Biophysical Journal, 46(1):15-27.

- 27 -

2.9

Using indexed variables in a linear array of cells

Features: • Indexed variables (arrays) • Plotting as a function of position (index) • Batch Runs Description of the model:

At time = 0, a linear array of n +2 cells (compartments) begin producing metabolite M at a constant rate q. M diffuses between cells with transfer coefficient k. The two end cells are slit open so that they are in equilibrium with an infinite bath where M = 0. Cell Concentrations and volumes are denoted by C and v respectively. The model computes C and M in all cells for all time >0. Follow these steps: 1. Run the model. A family of 52 kinetic curves corresponding to the 52 compartments is plotted. 2. Replot the results as a function of position, double click anywhere in the x-y plane of the graph (or select Choose Variables from the Graph menu.) The following dialog appears

Click on the X Axis popup and select the index [i] for the x axis. Click on OK. MAdonna now plots the concentration as a function of position in the array at the STOPTIME (presumably when the system has reached a steady state). Click OK and run the model 3. Plot a family of concentration-position profiles that show the transition from initial to steady state, click on Batch Runs in the Parameters menu and set it up as shown. In this case we want to make several automatic runs for different durations, so we select STOPTIME in the Parameter popup. We will display 10 runs with durations ranging from STOPTIME =0 (initial conditions) to 500 time units. The intermediate values are automatically calculated and displayed in the Values: scroll window. The intermediate values can be computed (Series Type:) on a linear (Arithmetic) or Geometric scale. However under Geometric, the initial value must be non-zero to make any sense. No averaging is required so the Keep Runs Separate radio button is activated.

- 28 -

4. Click OK. The concentration profiles are calculated and displayed. The time corresponding to each profile is displayed in the legend. (The legend is toggled by a L button in the top tool bar of the graph window.) 5. You can change n, the number of compartments ad lib by moving the slider. Model Equations Refer to the "Array Equations" section in the Equation Help window (Help menu) for an explanation of array syntax. The kinetic equations for the inner 50 compartments are defined in terms of the number of moles, M, in each compartment: d/dt (M [1..n]) = a + p*(C [i-1] -C[i] ) +p*(C [i+1] -C[i] ) Init M[1..n] =0 With equal volumes V in each compartment, concentrations are defined as C [1..n] = M[i] / V and the first and last compartments are reserved for the boundary conditions: C[0] = 0 C[n+1] = 0 -------------------------------------------------------------The tutorials are designed to introduce many of Berkeley Madonna's features by leading you step by step through a computation. In many instances the examples will be taken from our laboratory course, which can be accessed on the web at: http://mcb.berkeley.edu/courses/mcb137/.

- 29 -

2.10 Monte Carlo Model of a single Na channel channel under a patch clamp Feature: • Demonstrates Monte Carlo Simulation generated by using the AVERAGE feature of Batch Runs • Demonstates the use of difference equations using the NEXT statement Description of the model: A single sodium channel is observed using a patch clamp. A step in membrane voltage from Eo to E is applied, and currents through the channel are recorded for several milliseconds. Individual runs appear chaotic, but when many runs are averaged they reproduce the macroscopic currents predicted by the Hodgkin Huxley Equations. The model computes the probability of finding an open channel, which is proportional to the sodium current across the membrane. Follow these steps: 1. Run the model several times to see individual runs. Not much happens; occasionally you may see the channel open for a brief period. 2. Run the same experiment 1000 times and take the average of these runs: a. Click on the Batch Runs... entry in the Parameter menu. (see figure below) No variation in parameters is involved in this operation so that the Parameter: pull down is left showing "None". b. Set the # of Runs to 1000. c. Under Mode: activate the radio button entitled Compute Mean d. Click OK

3. Display the macroscopic Hodgkin Huxley prediction. Toggle the m3h_HH button on the graph. (Averaging the 10.000 simulation runs would give even better agreement.) Model Equations

Each channel is guarded by a gate which must open in order for the channel to pass ions. Each gate has 4 "subgates", 3 m gates, y1, y2, y3, and one h gate, y4. All subgates have to be open in order for - 30 -

the channel to be open. y is a variable giving the status of a particular subgate (i.e. y = 1 for open, y = 0 for closed). The variable m3hy = y1*y2*y3*y4 gives the probability of finding the channel in the open configuration. Behavior of each subgate is determined as follows: 1. An initial state, y, is calculated from qmo, the steady state fraction of open gates obtained from Hodgkin Huxley equations: INIT y = if RANDOM(0,1) > qmo then 1 else 0 where RANDOM(0,1) is a Madonna function that returns a random number between 0 and 1 2. An initial time, Tg, for the next event, is calculated from the assumption that events follow a Poisson distribution with probability am*dt for opening (and bm*dt for closing) within the interval, dt. The am and bm are the voltage dependent rate constants for opening and closing given by the Hodgkin Huxley equations. INIT Tg = -((1-y)/am + y/bm)*LOGN(RANDOM(0,1)) Note that: (1-y)/am + y/bm) = 1/am for y = 0 = 1/bm for y = 1. If p = the probability of the open event occurring at time t, then p = exp(-am*t), and t = (1/am)*LOGN(p). The usual assumption of equal apriori probabilities sets p = RANDOM(0,1) 3. Subsequent states as well as event times are generated in the same manner as soon as the next event occurs NEXT y = IF TIME = 0 {Value of "y_real" delayed by "delay_time":} y_delayed = delay(y_real, delay_time) flow a_real = y_real flow a_delayed = y_delayed init a_real = 0 init a_delayed = 0 delay_time = .5

- 33 -

3.3

Difference Equations

This model illustrates how you can compute the mean and standard deviation of a variable over time using difference equations. As the model runs, the "sum_x" and "sum_x2" variables are used to accumulate the sum and sum-of-squares of x. The mean and standard deviation are then computed from these variables. If you scroll to the end of the table, you can see that the NORMAL built-in does indeed generate a sequence of pseudo-random numbers with the specified mean (3) and deviation (0.8).

starttime = 0 stoptime = 10000 dt = 1 x = normal(3, 0.8) next sum_x = sum_x + x next sum_x2 = sum_x2 + x*x avg = if n > 0 then sum_x / n else 0 std = if n > 0 then sqrt((sum_x2 - sum_x * avg) / n) else 0 init sum_x = 0 init sum_x2 = 0 init n = 0 next n = n + 1

- 34 -

3.4

Floating Point Exceptions

This model illustrates how Madonna handles floating-point exceptions during model execution. In this case, the value of y overflows at times of -0.5 and +0.5 (approximately). Run the model and click "No" when the error message is displayed. You'll see that Madonna stopped the run at the time the overflow occurred. This is easier to see by switching to table view. Now go to the Model Settings dialog and uncheck the "Stop On Exception" box and run the model again. Note that even though an error occurred, Madonna continues the run until STOPTIME is reached. Also note how the plotter handles results that are out of range: it simply skips them.

- 35 -

Multinomial Function (Arrays) Implementation of MULTINOMIAL(p[*], N) This is like the BINOMIAL built-in but generates a bunch of outputs from a bunch of probabilities. Note that the probabilities don't need to add up to one. Instead, this algorithm scales them by their sum. This is a good demonstration of array indexing techniques.

N = 1000 {Number of trials} K = 5 {Number of probabilities} p[1] = 1 {List of probabilities} p[2] = 2 p[3] = 3 p[4] = 6 p[5] = 8 {The vectors "s" and "x" are used internally to compute partial sums, the result is placed in vector "y".} s[K] = p[K] s[K-1..1] = p[i] + s[i+1] x[0] = 0 x[1..K] = x[i-1] + binomial(p[i] / s[i], n - x[i-1]) y[1..K] = x[i] - x[i-1] {Sanity check: sum of y's elements should equal N} z = arraysum(y[*])

- 36 -

3.5

NETFLOW Built-in

This model illustrates why you need the NETFLOW built-in to compute the actual change in a reservoir when using higher-order integration methods. "flow_1" is the net flow into the reservoir "y" computed manually. Note that it matches the value of the NETFLOW built-in (stored in "flow_2") when using Euler's method, but it underestimates the flow for higher-order methods like Runge-Kutta 2.

STARTTIME = 0 STOPTIME = 3 DT = 0.25 INIT y = 1 FLOW y = dy dy = y ;This makes y increase exponentially init flow_1 = 0 next flow_1 = dy * DT flow_2 = netflow(y)

- 37 -

3.6

Phase Plane

This model illustrates the "Initial Conditions" (Ic) button in the graph window. By setting the X and Y axes to reservoirs with constant initializers, we can set their initial values simultaneously by clicking in the plot. Also, note that auto-scaling is turned off for both axes so that the scales don't change from run to run.

method rk2 dt = .01 init x = 0 init y = 1 flow x = -y flow y = x+k*y k = -0.3

- 38 -

3.7

Relations and Logical Operators

This model illustrates usage of Madonna's relational and logical operators. The resulting output looks like some kind of digital circuit. The PULSE function can be found on the Equation Help page: PULSE(v, f, r):

Impulses of volume v, first time f, repeat interval r

init d0 = 9 flow d0 = pulse(1, 1, 2) + pulse(-1, 2, 2) init d1 = 7 flow d1 = pulse(1, 2, 4) + pulse(-1, 4, 4) Tor = 5 + (d0 >= 9.5 or d1 >= 7.5) Tand = 3 + (d0 >= 9.5 and d1 >= 7.5) Tnor = 1 + not(Tor >= 5.5)

- 39 -

3.8

Stiff Systems

This model illustrates how the "Rosenbrock (stiff)" integration method is better suited to stiff problems. This example is taken from Numerical Recipes in C. Try running the model with the default Autostepsize algorithm. It takes over 51000 time steps. And if you try relaxing the tolerance (try 1e-8) you'll see the model output is starting to destabilize (zoom in on y3). Now switch to the Rosenbrock (STIFF) method and run it again. Same results with only 96 time steps! And you can reduce the tolerance and still get accurate solution with even fewer steps.

method auto starttime=0 stoptime=50 dt = 2.9e-4 { Initial minimum dt } dtmax=50 { allow huge steps } tolerance=1e-10 init df = 0 next df = dt / delta_t delta_t = dt init istep = 0 next istep = istep+1 f1 = -.013*y1 - 1000*y1*y3 f2 = -2500*y2*y3 flow y1=f1 flow y2=f2 flow y3=f1+f2 init y1=1 init y2=1 init y3=0

- 40 -

3.9

Vector Finals

This model illustrates the use of the "vector final" option in the plotter. It collects the results of the BINOMIAL built-in into bins and displays the final bin values in a single plot.

starttime = 0 stoptime = 10000 dt = 1 dtout = 1000 b = binomial(p, n) p = 0.5 n = 100 init a[0..n] = 0 next a[0..n] = a[i] + (b = i)

- 41 -

4 4.1

HOW DO I Control Plotted Variables

The model for the reaction sequence A → B → C → D → E → F → G → H → I → J → K → L → M → N is now active. • Run the model. Madonna automatically plots the first two variables in the equation window (in this case A and B). The first is scaled on the left axis, the second on the right. Buttons for the first eight variables (A...H) are placed on the left bottom margin of the display. • Click on the buttons to toggle the graphs of A,B, and C ... H. Toggle C on. • Use shift-click B to change the scaing axis for B from right to left axis. • The left axis is scaled to A because it has the largest value of any of the variables that refer to that axis. Use control-click C to change the left axis scaling from A to C. • For more general control of plotting variables double click anywhere in the x-y plane to bring up the following Choose Variables dialog (also available from the Graph menu).

All of the model variables are listed in the Variables list. In the illustration we have selected J which is not plotted (its numerical values were not saved during the last run). • Clicking the Add button (or double clicking the selection) will move J into the Y axes column and a J button will be added to the display. The next time the model is run, J will be available for plotting. • In the Y axes column we have selected F. Notice the check boxes below. The Right Axis is checked, indicating that F will be scaled on the Right. Removing the check will cause F to be scaled on the left. The Visible box is unchecked. This indicates that F will not be plotted until you toggle its button on the graph. • To remove F from the Y Axes list, click on Remove (or double click on F). Now the F button will disappear from the graph and F will be moved to the Variables list. • Note the pull down menu labelled "X Axis". The default variable plotted on the X axis is Time. However, All variables are listed in the pull down menu, and any of these can replace Time on the X axis.

- 42 -

METHOD RK4 STARTTIME = 0 STOPTIME = 20 DT = 0.02 d/dt (A) = - ka*A d/dt (B) =ka*A - kb*B d/dt (C) =ka*B - kb*C d/dt (D) =ka*C - kb*D d/dt (E) =ka*D - kb*E d/dt (F) =kas*E - kbs*F d/dt (G) =kas*F - kbs*G d/dt (H) =ka*G - kb*H d/dt (I) =ka*H - kb*I d/dt (J) =ka*I - kb*J d/dt (K) =ka*J - kb*K d/dt (L) =ka*K - kb*L d/dt (M) =ka*L - kb*M d/dt (N) = kb*M init A =100 init B = 0 init C = 0 init D = 0 init E = 0 init F = 0 init G = 0 init H = 0 init I = 0 init J = 0 init K = 0 init L = 0 init M = 0 init N = 0 ka =1 kb =1 kas=.1 kbs=.1

- 43 -

4.2

Define and Use Sliders

The model for simple chemical reaction A → B → C is now active. • Run the model • Open the Define Sliders dialog (in the Parameters menu). • Double click on ka and on kb to add them to the Slider list. • Slider behavior can be tailored to your specific needs by utilizing the entry fields on the right hand side of the dialog. Movement of the slider can produce either linear or logarithmic increments in the value of the parameter; the minimum and maximum parameter values corresponding to the two ends of the slider can be specified; and the sensitivity of slider movements can be adjusted by changing the increment. (For example, to constrain the parameter values to integers, set the increment =1) • Click OK. • Move the sliders back and forth: when you release the mouse, Madonna automatically runs the program with the current parameter values. The double arrow (just to the right of 10x) on the slider window operate just as in a scrolling window; they "fine tune" the slider movements. Clicking on them increases (or decreases) the parameter by one increment. • Each slider has a check box labeled 10x which reduces the range of the slider by a factor of ten around its current value. This feature makes it easier to zero in on a desired value. Note that this 10x mode is cancelled if you make any changes via the Define Sliders dialog or recompile your model. • Double click within the text area of the sliders to re-open the Sliders dialog where you can add sliders, adjust the current slider limits etc.

METHOD RK4 STARTTIME = 0 STOPTIME = 10 DT = 0.02 d/dt (A) = - ka*A d/dt (B) =ka*A - kb*B d/dt (C) = kb*B init A =100 init B = 0 init C = 0 ka =1 kb =1

- 44 -

4.3

Display Numerical Values

The model for the reaction sequence A → B → C is now active. • Run the model. • Retrieve data directly from the plot: Activate readout mode by depressing the Readout button, then drag the mouse around inside the data area of the plot. As you drag, the "crosshair" cursor tracks the mouse and the corresponding values are displayed in the information area at the top of the window. (When Fast Fourier Transform mode is active, the readout displays both frequency and period corresponding to the X axis coordinate.) Note that you must deactivate readout mode in order to use the "zoom in" feature. • Toggle the plot into numerical tables: Click on in the graph toolbar to display the contents of the graph in tabular form. (For superimposed plots, note that the TIME column reflects the time for the run shown in the rightmost column of the table. If there are other runs with different time scales, their values will be mapped to the times shown in the TIME column using linear interpolation.) Variables can be added and deleted from the table in the same way that they are added and removed from the plots (Use the buttons and/or double click on the contents of the table - see "How do I... Control Plotted Variables" in the Help menu.) To return to the plot, simply click on again. • Export Numerical Data Tables can be saved as tab-delimited text or comma-delimited text (CSV) files by activating the window containing the table and choosing Save Table As from the File menu. Tables can also be copied to the clipboard using the Copy Table command in the Edit menu. Tables are always copied in tab-delimited text format.

METHOD RK4 STARTTIME = 0 STOPTIME = 10 DT = 0.02 d/dt (A) = - ka*A d/dt (B) =ka*A - kb*B d/dt (C) =kb*B init A =100 init B =0 init C =0 ka =1 kb =1

- 45 -

4.4

Make a Phase Plane Plot

The model for a non-linear oscillator is now active. • Run the model • Double click anywhere in the x-y plane of the graph (or select Choose Variables from the Graph menu.) The following dialog appears...

• Remove Z from the Y axis. • Replace TIME on the x axis with Z by selecting Z in the drop down X Axis menu. • Click OK. You will see a phase plane plot starting from the initial condition and moving into a closed loop limit cycle. Change initial conditions. • Toggle the (overlay) button on the graph tool bar. (This allows results to be compared, but it is not required.) • Toggle the (initial conditions) button. The graph cursor changes from to . • Place the cursor anywhere in the Y-Z plane and click the mouse. Madonna interprets this as a new set of initial conditions, for Y and Z, given by the coordinates of the point that you have selected. Madonna automatically runs the program with these new conditions; numerical values are displayed in the graph window, just above the graph. • Repeat this at several different positions. In each case you will see the plot converge on the same limit cycle

- 46 -

4.5

Perform a Curve Fit

A. Example: Single Curve (see part C below for general description) The model for simple chemical reaction A → B → C is now active. • Run the model. Curve fitting requires an initial guess of all unknown parameters. In our case we guess ka = .1 and kb = .3. These are entered in the model as though they were known (see Equation Window). Running the model shows the results of our guess. • Import data a. Select Import Dataset from the File menu. The standard Open dialog will appear b. Open the enclosed data file entitled "Bdata" (contained in the "How do I" folder within the Berkeley Madonna folder). This is a simple 2 column text file (Data can be text or spread sheet files with an unlimited number of columns). The dialog shown below will appear. It describes the data file, allows you to rename it, and allows you to choose which column represents the independent (x axis) variable and which represents y. Numerical data files in Berkeley Madonna are referred to by a # prefix. The default name is shown below. c. Click OK.

The data is now super imposed on the plot. Notice the large discrepancy between the data and the model run with the guessed parameters. • Curve fit a. Click Curve fit in the Parameters menu. The dialog shown below appears. b. Double click on the unknown parameters ka and kb in the Available: list to move them over into the unknown Parameters list. c. Use the Fit Variable pull down menu to select B. d. Check that #Bdata is entered in To Dataset . e. Click OK. You are finished! Parameter values used in the plotted curve fit can be found in the current Parameter window. (See below and, even better, the Users manual for more options on curve fitting. For example, the right hand side of the dialog can be used to limit the permissible parameter values, multiple curves can be fit simultaneously, rms error can be retrieved, etc.)

- 47 -

B. Example: Multiple Curves Berkeley Madonna will also fit data to multiple curves; in this example of A ∅ B ∅ C, we import data for A, B, and C and look for parameters to fit all three curves simultaneously. • Import data a. Select Import Dataset from the File menu. b. Open the enclosed data files entitled "Adata", "Bdata" (if it has not already been imported) and "Cdata" (all contained in the "How do I" folder within the Berkeley Madonna folder). c. Click OK • Curve fit a. Click Curve fit in the Parameters menu. The dialog shown below appears. b. Double click on the unknown parameters ka and kb in the Available: list to move them over into the unknown Parameters list. c. Click on the Multiple Fits checkbox. d. Use the Fit Variable pull down menu to select A. e. Select #Adata from the To Dataset list. f. Click on the Add button to list the dataset in the multiple fits list g. Repeat steps d, e, and f, for B and C h. Click OK. You are finished! Note the values of the fitted parameters ka, and kb. The correct values (i.e. those that were used in the generation of the "empirical" data) are ka = 1.00, kb = 1.00.

- 48 -

C. General Description (from User's Guide p 45) Curve Fitter Berkeley Madonna can automatically find the values of one or more parameters in your model that minimize the deviation between your model's output and a dataset. Before you use this feature, you should be able to run your model and be confident that it is working properly. You should already have determined integration parameters that give accurate results while minimizing execution time. To perform a curve fit, choose Curve Fit from the Parameters menu. You will see a fairly involved dialog box. Here's what you do with it: • Specify which variable in your model you are trying to fit to the external data. • Specify the dataset you want to fit the variable to. If you haven't imported the dataset yet, you can do so by clicking the Import Dataset button. See Importing Datasets on page 15 for details. • Specify one or more parameters in your model which you want to solve for. To do this, choose a parameter from the available list on the left and click Add. If you change your mind, click Remove. • For each parameter you have chosen, you must provide two initial guesses. To set guesses for a parameter, select its name in the parameter list and edit the values displayed in the Guess #1 and Guess #2 boxes. • You must also provide a minimum and maximum value for each parameter. Both guesses must be within this range. The curve fitter will not set the parameter outside of this range. • Specify the fractional tolerance you desire in the solution. For example, if you want the parameter(s) solved to three significant figures, set the tolerance to 0.001. The tolerance must be greater than zero. • Once you are satisfied with your setup, click OK to begin the curve fit. Berkeley Madonna will run your model repeatedly until it finds a solution. If you get tired of waiting, click the Stop button to abort the fit. • When the fit is complete, Berkeley Madonna leaves the parameters set to the values that give the best fit. It then runs your model and plots the fit variable and dataset in a graph window. While the curve fit is in progress, Berkeley Madonna displays the RMS deviation between the data and best run so far. The deviation is root mean square of the differences between individual data points in the dataset and the corresponding points in the run.

- 49 -

When the curve fit operation complete, Berkeley Madonna immediately closes the Running dialog. If you want the dialog to remain open after the curve fit completes, check the Pause After Curve Fit option in the General page of the Preferences dialog. This enables you to observe the RMS deviation of the best fit. To close the dialog, click the Done button. Multiple Fits By default, the curve fitter fits a single variable in your model to a single imported dataset. However, it can simultaneously fit multiple variables to multiple datasets. To do this, check the Multiple Fits box in the Curve Fit dialog and add variable-dataset pairs to the list of fits using the Add button. Each variable-dataset pair has a relative weight which allows you to control the importance of each pair. The default weight is one, but it can be changed by selecting the pair and entering a new weight in the box below the list of fits. Berkeley Madonna multiplies the deviation of each pair by its weight, then sums these products to compute the overall deviation to minimize. Note that when the Multiple Fits box is checked, the selections in the Fit Variable and To Dataset controls have no effect on the curve fit operation; they are used simply to define variable-dataset pairs to add to the Multiple Fits list. Specifying Initial Guesses There are subtleties involved in specifying initial guesses for each parameter. First, the two guesses must be different. If one or more parameter's initial guesses do not meet this criterion, the curve fit will stop before the best solution is found. When you add a parameter to the list, Berkeley Madonna proposes default guesses. Guess #1 is set to 0.5 times the parameter's value. Guess #2 is set to 1.5 times the parameter's value. However, if the parameter's value is zero, the guesses are set to -1 and +1. Floating-point Exceptions Some models are very touchy when you vary their parameters. That is, they will not run at all! When this happens, you get a floating-point exception. As the curve fitter runs, it may change the value of a parameter in such a way that the model won't run. The curve fitter can handle this case by backing off, but when it first starts up it must be able to run your model with any combination of your initial guesses. If your model won't run with one of these combinations, the curve fitter won't be able to start. You'll know when this happens if you get an error while the curve fitter is in the initialization phase.

- 50 -

4.6

Plot a Family of Curves

The model for simple chemical reaction A → B → C is now active. • Run the model • Click on Batch Runs in the Parameters menu. • Set up the dialog as shown below. In this case we want to make several automatic runs for different values of ka, so we select ka in the Parameter pull down menu. We will display 10 runs with values ranging from ka = 0.1 to 1.0. The intermediate values are automatically calculated and displayed in the Values: scroll window. With Arithmetic checked, the intermediate values are separated by equal increments. With Geometric checked, the values form a geometric sequence.

• Click OK. The concentration profiles are calculated and displayed. • Toggle the L button. The ka corresponding to each profile is displayed in the legend. • Toggle A off the graph so that you can see the unencumbered behavior of B. The same can be done for A and C.

METHOD RK4 STARTTIME = 0 STOPTIME = 10 DT = 0.02 d/dt (A) = - ka*A d/dt (B) =ka*A - kb*B d/dt (C) = kb*B init A =100 init B = 0 init C = 0 ka =1 kb =1

- 51 -

4.7

Solve Boundary Value Problems

A. Example (see part B below for for general description) • Example: Steady state radial diffusion from a cylindrical surface with radius = 5 µ outwards into an infinite sink located at r=100 µ. The steady state diffusion equation in cylindrical coordinates is:

The constant source at the surface and the infinite sink at 100 µ provide the following boundary values: C(5) = 1 C(100) = 0 • The independent variable is now r, rather than the default TIME. This change is inserted in the program by the statement: RENAME TIME = r

(see Equation window)

• Higher order derivatives can be processed by Madonna if they are written with "prime" notation, and if the highest order derivative is isolated and placed on the left hand side of the equation. In "prime" notation the diffusion equation becomes D(C'' +C'/r) =0, and this is entered as: C'' = - C'/r • The initial condition for C' is not known. It is arbitrarily set at -2. The BVP solver will look for the appropriate initial condition; it will keep changing this condition and running the model until it finds one where C = 0 at r = 100. This is activated by the Boundary Value ODE dialog (Model menu - Modules popup). The setup in the boundary value dialog box is shown below. To set the boundary condition: 1. Select C from the Set menu in the drop-down menu at the top of the dialog 2. Fill in the 2 edit boxes that lie immediately below ("=" and "at X") 3. Click on the Add button. This will automatically enter the condition in the Boundary Conditions list. 4. At the bottom of the dialog, choose which parameter to vary in order to satisfy the boundary condition. In our case, we select INIT C' since we have assigned an arbitrary value to it. 5.Click on the Add button. This will automatically enter the parameter in the Unknowns list. 6. The Minimum and Maximum edit boxes can be adjusted to narrow the search by setting upper and lower limits to the unknown.

- 52 -

B. General Description (from User's Guide p 47) Boundary Value Solver Berkeley Madonna can solve boundary value problems by varying one or more parameters such that specified boundary conditions are met. The requirements for using the boundary value solver are similar to those for the curve fitter. However, instead of using imported datasets, the boundary value solver uses one or more boundary conditions that you specify in the Boundary Value ODE dialog box. To solve a boundary value problem, choose Boundary Value ODE from the Modules submenu in the Model menu. Then, complete the dialog box as follows: • Specify a boundary condition by: (a) choosing a variable from the Set control, (b) entering the target value in the = field, (c) entering the time at which the target value should be reached in the at X = field. Then, click the Add button. The boundary condition will be displayed in the form var(x)=y in the Boundary Conditions list on the right. Repeat this procedure for each boundary condition. • To remove a condition you don't want, select it from the list and click Remove. To change a condition, remove it, edit the appropriate fields, then add it back again. • Specify a parameter to solve for by choosing one from the Available Parameters list and clicking the Add button. The parameter appears in the Unknowns list on the right. • Adjust the minimum and maximum values for the parameter so that they are just wide enough to cover the range where you expect the solution to be found. Berkeley Madonna will not search for solutions outside of these limits. If the limits are too large, it will take much longer to find a solution or may not be able to find one at all. • Repeat steps 3 and 4 for each parameter you want to solve for. Note that the number of parameters must be the same as the number of boundary conditions specified in step 1. • To remove a parameter you don't want, select it from the Unknowns list and click Remove. • Specify the relative accuracy you desire in the Tolerance field. The default is 0.001 (one part per thousand). Once everything is set up, start the solver by clicking OK. The solver will run until it finds a solution that meets all of the specified boundary conditions. If a solution cannot be found, the solver will run indefinitely. You can stop it at any time by clicking the Stop button. When the solver finds a solution, it leaves the parameters set to the values that satisfy the boundary conditions and runs your model, plotting the variables used in the boundary conditions. If you stop the solver before it finds a solution, it leaves the parameters set to the values that came closest to the solution and does not run your model. Sometimes the solver will complete as if it had found a solution, but when you examine the results of the run it appears that the boundary conditions you specified aren't met. This happens because the solver doesn't actually require all conditions be met exactly (if it did, it would never finish due to the limited precision of floating-point numbers on a computer). Instead, the solver considers its work done when the relative adjustments it makes to each parameter are less than the specified tolerance. If you model doesn't respond smoothly to small parameter changes close to the solution, the solver may think it's close to a solution when in fact it is not. To avoid this problem, try using a lower tolerance in the Boundary Value ODE dialog. Or try using a smaller stepsize to improve the stability of your model's output.

- 53 -

METHOD RK4 STARTTIME = 4 STOPTIME = 100 DT = 0.02 RENAME TIME= r C'' = -c'/r init C = 100 init C'=-2

- 54 -

4.8

Chemical Kinetics Module

Transient State Michaelis Menten Kinetics Features: The chemical reactions module is a template that can be appended to Berkeley Madonna's models. It allows you to enter reactions in ordinary chemical notation; Madonna will automatically set up the corresponding kinetic equations based on the law of mass action and enter them into the equation window. The equation window is fully accesible for writing additional equations. Description of the model: Non steady-state Michaelis Menten model for enzyme kinetics is set up and integrated. Follow these steps: 1. Open the Chemical Reactions dialog by selecting it from the Model menu: Modules popup. 2. Using ordinary chemical notation, enter the first reaction in the Reactants: and Products: boxes as shown below. 3. Enter Kf, the forward rate constant, and Kr, the reverse rate constant, as shown. 4. Click Add, and the reaction will be recorded in the Reactions: list. 5. Enter the Initial Concentrations as shown.

- 55 -

6. Repeat the same operation for the second reaction:

7. Click OK. The following will appear in the Equation window:

The gray areas display Madonna's interpretation of the Chemical Reaction dialog. Change the STOPTIME to 300. (All parameters can be changed in the usual way with the Parameters window.) Click the Run button to see the results. It is instructive to plot the results against log (time). Double click on the x axis to open a dialog for making this change Important Notes: Contents of the gray areas can be supressed by the Show check boxes on the lower right in the Chemical Reactions dialog. Try them. Edit the gray areas. Double click on the gray areas and re-open the reactions dialog. Select the reaction to be edited from the reactions list and change it in the top Reactants: and Products: boxes. Once a reaction is changed it is entered with the Modify button. If the Add button is used Madonna will interpret it as an additional reaction. To remove a selected reaction use the Remove button.

- 56 -

To change the mass action interpretation of any reaction rate simply rewrite the rate in the Equation window. E.g. if the desired first reaction rate is given by f(S,E), then write RXN1 = f(S,E) in the Equation window, and it will supercede the reaction written in the gray area, provided that the "new" statement is entered after the "old" one. Overwriting the chemical reaction module In general, equations entered into the equation window following the gray area will take precedence. However, a differential equation cannot be replaced by an algebraic expression and vice versa.

- 57 -

4.9

Using Conveyors

Introduction A conveyor is defined like this: name = conveyor(input, tt) where name is the name of the conveyor, input represents material going into the conveyor and tt is the transit time. Each time the model is stepped, the input is added to the conveyor and its transit time is noted. Later, when the input's transit time has expired, it flows out of the conveyor. The value of the conveyor variable is the sum of all inputs in the conveyor multiplied by DT. The conveyor's outflow can be determined using the OUTFLOW built-in function: output = outflow(name) You can prevent the total amount of material in a conveyor from exceeding a specified capacity. To do this, you provide a capacity argument when defining the conveyor: name = conveyor(input, tt, capacity) Then, use the MAXINFLOW built-in to determine how much the conveyor can accept during a particular time step: input2 = min(input, maxinflow(name)) name = conveyor(input2, tt, capacity) Here, input2 is constrained to prevent the conveyor's capacity from being exceeded. Example A capacity-constrained conveyor named "C" is shown in the Equation window. Run the model. Notice that the first three input pulses are added to the conveyor, but the fourth input pulse (at TIME=10) is rejected because the conveyor has been filled to capacity. The next pulse (at TIME=13) is accepted because the very first pulse exited the conveyor at TIME=11 making room for another input. Try increasing the capacity of the conveyor to 4 or more. Now the conveyor can accept every input pulse.

METHOD Euler STARTTIME = 0 STOPTIME = 25 DT = 1 capacity = 3 tt = 10 input = pulse(1, 1, 3) input2 = min(input, maxinflow(C)) C = conveyor(input2, tt, capacity) output = outflow(C)

- 58 -

4.10 Use Indexed Variables (Arrays) Choose "Equation Help" from the Help menu and look over the section titled "Array Equations". Note 1: The symbols i, j, and k are reserved for the running index; they are the only symbols that can be used for this purpose and are available only within the right-hand side of an array equation. Note 2: The left-hand side of an array equation uses the range operator (two dots "..") to define the extent of each dimension (the first and last index). Example: • The model for the reaction sequence A1 ∅ A2 ∅ ... ∅ Ai ∅ ... ∅ An is now active. It is written in array notation. The first reactant A1 is not produced and the last reaction is not consumed so they are handled separately. The rest of the reactions i = 2...(n-1) are covered by the equation d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i] where k[i] = the rate constant for the dissipation of the ith reactant • The first and last reactions are d/dt (A[1]) = -k[1]*A[1] d/dt (A[n]) = k[n-1]*A[n-1] • The redundant initial conditions also use indices, as do the rate constants (see Equation window) • Run the model. The entire indexed set (designated by A[ ] is plotted). Use the legend to distinguish the different plots • Notice the last two statements which define the variable Ax = A[m] with m arbitraily set to 6. Clicking the Ax button on the graph will plot A[m] = A[6] for the default. We then define a slider for m to make it easy to obtain individual plots for any member of the set.

METHOD RK4 STARTTIME = 0 STOPTIME = 20 DT = 0.02 d/dt (A[1]) = -k[1]*A[1] d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i] d/dt (A[n]) = k[n-1]*A[n-1] init A[1] = 100 init A[2..n]=0 k[1..n] =1 n = 14 Ax = A[m] m= 6

- 59 -

4.11 Using Multidimensional Arrays Read the "Array Equations" section in the Equation Help window (Help menu) for basic information about array equations. This model illustrates flow of material over a surface. The surface "S" is modelled as a 10x10 array of cells. Each cell is a differential equation that accepts material from three of its neighboring cells (above, left, and above-left diagonal). The model begins with material in the top-left cell only. As time advances, material flows in a southeast direction over the surface and eventually falls off the right and bottom edges. The last equation in the model copies the row of cells at the bottom of the sheet into a onedimensional array (vector). This slicing technique lets you plot only the portion of the array you're interested in. Note that the size of the matrix can be easily changed by adjusting the values of M and N in the Parameter window.

{dimensions of matrix} m = 10 n = 10 {initial material in top-left corner} init S[0..m,0..n] = 0 ;Sets entire matrix to zero init S[0,0] = 10 ;Changes top-left cell {propagation constants} a = 2/(sqrt(2)+4) ;Horizontal/Vertical b = 1/(sqrt(8)+1) ;Diagonal {flow in southeast direction only} d/dt(S[1..m,1..n]) = a*S[i,j-1] + a*S[i-1,j]+b*S[i-1,j-1] - S[i,j] d/dt(S[0,1..n]) = a*S[i,j-1] - S[i,j] d/dt(S[1..m,0]) = a*S[i-1,j] - S[i,j] d/dt(S[0,0]) = -S[i,j] {cross section across bottom of sheet} slice[0..n] = S[m,i]

- 60 -

4.12 Using Ovens Introduction An oven is defined like this: name = oven(input, ct, ft, cap) where name is the name of the oven, input represents material going into the oven, ct is the cook time, ft is the fill time, and cap is the capacity. The capacity argument is optional; if omitted, the oven has unlimited capacity. Ovens accumulate input (filling state), hold these inputs for a specified time (cooking state), and release them (emptying state). The value of the oven variable reflects the amount of material in the oven multiplied by DT. An oven starts out in the idle state. When its input becomes nonzero, it begins filling. While the oven is filling, the MAXINFLOW built-in indicates how much material the oven can accept without exceeding its capacity limit: max_input = maxinflow(name) The oven switches to the cooking state when one of two things happen: the fill time (sampled when the oven switches to the filling state) has expired or the oven's capacity has been reached. While the oven is cooking, it does not accept input and its value remains constant. The oven enters the emptying phase when its cook time has expired. The emptying phase lasts for only one time step. During this time, the entire contents of the oven are released. This can be observed using the OUTFLOW built-in function: output = outflow(name) The oven begins filling again during the emptying phase if there is a nonzero input present during this time. If no input is present, the oven switches back to the idle state and waits for input. The state of an oven can be determined using the OSTATE built-in: state = ostate(name) State values are 0 (idle), 1 (filling), 2 (cooking), and 3 (emptying) Example A model containing a reservoir "R" and an oven "O" is shown in the Equation window. The oven can accept up to half of the material from the reservoir during any given time step. Run the model and examine the output. Look for the following events: • TIME=0: The oven is idle and waiting for input. The reservoir has a value of 10 so the oven can accept up to 5. Since the oven's capacity is only 4, it fills up during this time step. • TIME=1: The oven begins cooking because it is full. • TIME=4: The oven finishes cooking and empties its contents. It also begins filling because an input is available. • TIME=5: The oven continues filling and reaches its capacity. • TIME=6: The oven begins cooking its second batch. • TIME=9: The oven finishes cooking and empties. It also begins filling again because an input is available. • TIME=10-12: The oven continues to fill. - 61 -

• TIME=13: The oven begins cooking its third batch. This isn't a full batch because the oven's fill time expired before its capacity limit was reached. • TIME=16: The oven finishes cooking the third batch and empties it.

METHOD Euler STARTTIME = 0 STOPTIME = 20 DT = 1 capacity = 4 init R = 10 d/dt(R) = -input input = min(R/2, maxinflow(O)) O = oven(input, 3, 4, capacity) output = outflow(O) state = ostate(O)

- 62 -

4.13 Using Queues Introduction A queue is defined like this: name = queue(input1, input2, ..., inputN) where name is the name of the queue and input1 through inputN represent items entering the queue. During each time step, a queue examines each of its inputs for nonzero values. For each nonzero input, an item with this value is placed in the queue. The value of the queue variable is the sum of all items on the queue multiplied by DT. Items are removed from the queue using the QPOP and QPOPS built-in functions: output = qpop(name, maxel, maxamt) output = qpops(name, maxel, maxamt) These built-ins remove ("pop") items from the queue in first-in first-out order. Up to maxel items are removed such that their sum doesn't exceed maxamt. The QPOP function removes only whole items from the queue whereas the QPOPS function can remove part of an item in order to reach the maxamt limit. Inputs are placed on the queue before the QPOP/QPOPS built-ins are evaluated which means they can be removed during the same time step that they are added. In this case, the queue's value is not affected since the items don't spend any time waiting in the queue. Example A model containing a queue "Q" feeding two ovens "X" and "Y" is shown in the Equation window. If you are not familiar with ovens, examine the "How Do I Use Ovens" model before proceeding. The queue has two pulse trains as inputs. When oven X is filling, it takes one whole item from the queue. When oven Y is filling, it takes as much as it can from the queue, splitting an item if need to reach the oven's capacity. Since the QPOP expression for oven X is listed before the QPOPS expression for oven Y, oven X gets first crack at the queue's contents. Run the model and note these events in the output: • TIME=1: A single input is added to the queue. Oven X immediately removes it and begins filling. • TIME=2: Two inputs are added to the queue. Oven X is still filling, so it removes the first element. Since there is another item in the queue, oven Y removes it and begins filling. • TIME=3: A single input pulse is added to the queue. Since both ovens have stopped filling (their fill times have expired), this input remains on the queue as reflected in the queue's value during the next time step. • TIME=4: Two more inputs are added to the queue and nothing is removed, so the queue now contains the following items from oldest to newest: 1, 1, 2. • TIME=5: Another input is added to the queue. Oven X has finished cooking and is now able to fill. Therefore, it withdraws the oldest item in the queue (1) and begins filling. Notice that the queue's value hasn't changed since the item removed is balanced by the item added during this step. The queue's contents are now (oldest to newest): 1, 2, 1. • TIME=6: Two more inputs are added to the queue. Oven X is still filling, so it removes another item from the queue. Oven Y is finished cooking and begins filling by withdrawing as much as it can from the queue. Since it has a capacity of 3.5, it removes the two oldest elements (2 and 1) and half of the third element (0.5). If the QPOP built-in had been used instead for oven Y's input, the third element would not have been split (try it). The queue's contents are now: 0.5, 2. - 63 -

• TIME=7: A single input is added to the queue. The ovens have started cooking again, so no items are removed from the queue. The queue's contents are now: 0.5, 2, 1. Notice how the amount of material on the queue increases over time since the ovens aren't capable of processing the inputs as quickly as they are created. If oven X were able to accept up to two items from the queue instead of one, the queue wouldn't fill up over time (try it).

METHOD Euler STARTTIME = 0 STOPTIME = 20 DT = 1 in1 = pulse(1, 1, 1) in2 = pulse(2, 2, 2) Q = queue(in1, in2) inX = qpop(Q, 1, maxinflow(X)) inY = qpops(Q, inf, maxinflow(Y)) X = oven(inX, 2, 2) Y = oven(inY, 3, 1, 3.5)

- 64 -

4.14 Using the Flow Chart Before starting, verify that the flowchart editor is functioning properly. To do this, choose About Berkeley Madonna from the Help menu and look for text that reads "Flowchart version XXX". If it says "Java not loaded", click the Load Java button. If Java fails to load, you need to install (or perhaps reinstall) Java support on your system. See the Read Me document for details. 1. Choose New Flowchart from the File menu. An empty flowchart window appears. 2. Position the mouse over the reservoir tool , press the mouse button, drag the reservoir to the flowchart, and release the mouse button. You've now placed a reservoir on the flowchart with the default name "R1". Note that the reservoir is colored red which means it's selected. 3. Change the name of the reservoir to "Population" by typing the new name. It is not necessary to first click the reservoir since it is already selected. 4. Click the flow tool , move the mouse over the Population reservoir, press the mouse button, drag the mouse to the right a few inches, and release the mouse button. This places a flow icon on the flowchart which drains the Population reservoir into an infinite sink. 5. Make sure the reservoir is connected to the flow. If the flow looks like this:

it's source end is not connected to the reservoir. In this case, drag the infinite source to the reservoir and drop it. This connects the flow's input to the reservoir instead of an infinite source. 6. You may find that the icons are not positioned the way you want. They can be repositioned by clicking and dragging. Note that Berkeley Madonna maintains the connection between the reservoir, flow, and sink. 7. Change the name of the flow to "Deaths" by clicking the flow if it is not already selected and typing the new name. on the flowchart approximately halfway between the reservoir and flow 8. Place a formula icon icons (use the same click and drag technique as you did for the reservoir). Change the name of the formula icon to "Death Rate". 9. Click the arc tool and create an arc going from the reservoir icon to the flow icon (use the same technique as you did for the flow). This creates a dependency relationship between the flow and the reservoir: the number of deaths depends on the population. 10. Create another arc going from the formula icon to the flow icon. This makes the number of deaths depend upon the death rate constant in addition to the population. 11. Double-click the Population reservoir. This opens a dialog box for the reservoir. In this dialog, you specify the reservoir's initial value. In this case, enter 1000 and click OK. Note that the question mark inside the reservoir icon disappears after the dialog has been dismissed. This means that the equations for this reservoir are fully specified. 12. Double-click the Death Rate formula to open it's icon dialog. Enter 0.02 for the formula's right-hand side and click OK. 13. Double-click the Deaths flow. Note that "Population" and "Death Rate" are required inputs. This means that these variables must be used in the right-hand side of the equation defining the icon's value. Enter the following formula for the right-hand side: "Population * Death_Rate". Note that instead of typing names of required inputs, you can double-click on the input's name in the list to insert that name. And you can use the calculator keypad to insert digits and common arithmetic operators. Click OK when done.

- 65 -

14. Your model's flowchart is now complete and should now appear something like this:

Now run your model by choosing Run from the Compute menu. Berkeley Madonna runs your model and displays the results in a graph window. 15. You can change the curvature of the dependency arcs by dragging their control points. To do this, first select an arc by clicking somewhere along its trajectory. The arc will turn red and a small square handle (control point) will appear. Click and drag the control point to change the arc's trajectory. 16. Icons and arcs can be removed by selecting them and choosing Clear (Macintosh) or Delete (Windows) from the Edit menu. You can also remove selected objects by holding down the control key while pressing the delete or backspace key. Try deleting the arc from the Population reservoir to the Deaths flow. Note that a question mark reappears in the flow icon. That's because the icon's equation refers to the Population variable which is no longer a dependent. You can easily reestablish this dependency by drawing a new arc from the reservoir to the flow. 17. Now that you've seen how easy it is to create a model, you should take a look at the equations Berkeley Madonna generated for you. To do this, choose Equations from the Model menu. Position the equations along side the flowchart so that you can see both windows simultaneously. 18. Select an icon in the flowchart with the mouse. When you select a single icon, Berkeley Madonna highlights its equation(s) in the equation window. This makes it easy to see an icon's equations without opening its icon dialog. 19. Berkeley Madonna can show you which icon corresponds to an equation in the equation window. To do this, position the equation window's insertion point (caret) within one of the equations, then choose Show Icon from the Edit menu. Berkeley Madonna activates the flowchart window and selects the corresponding icon. 20. Now try editing your equations in the equation window like you did in the first section of this tutorial. Berkeley Madonna won't let you. When working with a visual model (i.e., one that has a flowchart), you cannot edit the equations directly in the equation window. Instead, you edit each icon's equations using its icon dialog. The icon dialog can be opened by double-clicking an icon or choosing Icon Info.. from the Flowchart menu. You can also open the icon dialog directly from the equation window by double-clicking an equation. For example, double-clicking anywhere within the line "Deaths = Population * Death_Rate" opens the icon dialog for the Deaths flow. 21. Select the Population reservoir, type "Molecules" and press return. The name of the icon changes and Berkeley Madonna updates the name in all equations that depend on it. To see this in action, change the name again and keep your eye on the equation window as you press the return key.

- 66 -

• Run. Runs the model once. Same as choosing Run from the Model or Compute menus, except for parameter plot windows where this button performs a parameter plot run. • Lock. Prevents the window's contents from being changed when running the model or discarding runs in other windows. • Overlay Plots. When enabled, new runs are added to existing runs in the window. Otherwise, existing runs are discarded before new runs are stored. • Table. Displays the contents of the window in tabular form. Note that the TIME column reflects the time for the run shown in the rightmost column of the table. If there are other runs with different time scales, their values will be mapped to the times shown in the TIME column using linear interpolation. • FFT. Displays the contents of the graph window in the frequency domain. See Fast-Fourier Transform on page 15 of the user's guide. • Legend. Displays the legend. The legend shows the variable name and run number for each curve. Additional information in parentheses is included for Batch Runs: the parameter value for separate runs and labels for mean and mean±sd runs. You can place the legend anywhere in the graph window by dragging it with the mouse. • Parameters. Displays the parameter list. This list shows the value of each parameter when the model ran. If there is more than one run in memory and a parameter is not the same in all runs, the minimum and maximum values of the parameter are shown. You can position the parameter list with the mouse just like the legend. • Oscilloscope Mode. Displays the data like an oscilloscope. This button is present only if the "trigger" symbol is defined in your model. See Oscilloscope Mode on page14 of the user's guide. • Colors. Displays each curve in a different color. • Dashed Lines. Displays each curve in a different line style (solid or dashed). • Data Points. Draws a solid dot centered over each point in your data. If you have many closelyspaced data points, this has the effect of making the curves very thick. They will also redraw more slowly when this option is enabled. Tip: if you hold the control key while clicking this button, the lines connecting the data points will be hidden. • Grid. Shows the grid lines. • Readout. Displays X and Y values as you drag the mouse over the plot, instead of zooming in. See Readout on page 14 of the user's guide. • Initial Conditions. Enables you to set initial conditions by dragging the mouse over the plot. See Initial Conditions on page 15 of the user's guide.

- 67 -

• Zoom Out. Undoes the effect of the last zoom-in operation. If no zoom-in operations have been performed, this button is grayed. Note: Some of the graph window options can be turned on or off by default for new graph windows by editing the Graph Windows page of the Preferences dialog.

- 68 -