Author: Camron Fletcher

Equation Help, Examples & Guidelines

-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

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 -

- 49 -

- 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 -

- 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 -

• 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 -

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: