## Pn1 pn At [bi2F12pnF12pniB12pn1 1252

Euler's method is called an explicit method because p 1 can be written explicitly in terms of pj. Each time we use the above technique to update pj, we introduce a small error due to the finite size of At and round-off error. In the absence of round-off error, we can achieve any desired accuracy by decreasing At. However, there are some problems with this approach. Usually, the biggest problem is the amount of computer time required when we choose a very small At. However, the round-off error...

## Bard Ermentrout and Joel Keizer

Nonlinear ordinary differential equations are notoriously difficult or impossible to solve analytically. On the other hand, the solution to linear equations like those encountered in the kinetic model of the GLUT transporter can be expressed in terms of simple functions, and their behavior analyzed using standard results from linear algebra. In the first two sections, we cover some basic ideas in linear algebra and a review of power series. These are concepts that are needed in the sections...

## B110 Command Summary The Basics

Computes a trajectory with the initial conditions specified in the Initial Data Window ( T G ). computes a trajectory with the initial conditions specified by m(I)ce lets you specify many initial conditions ( Y lets you define a new 2D view ( V 11 2 ). Graphic stuff Postscript allows you to create a PostScript file of the current (F)it fits the window to include the entire solution ( W 11 F ). saves the state of XPPAUT( restores the state of XPPAUT from a saved . set file (

## Spatial Modeling

All of the models considered in previous chapters have relied on the implicit assumption that chemical concentrations are uniform in space. This assumption is reasonable when the region of space in which the reaction takes place is confined and quite small. However, there are many situations in which chemical concentrations are not uniform in space. A well-known example in which nonuniform distributions are crucial is the propagation of an action potential along the axon of a nerve fiber...

## B14 Solving the Equations Graphing and Plotting

Here, we will solve the ODEs, use the mouse to select different initial conditions, save plots of various types, and create files for printing. In the Main Window you should see a box with axis numbers. The title in the window should say X vs T, which tells you that the variable X is along the vertical axis and T along the horizontal. The plotting range is from 0 to 20 along the horizontal and 1 to 1 along the vertical axis. When a solution is computed, this view will be shown. Click on Init...

## Conservation Law in Multiple Dimensions

Consider a chemical species C whose concentration c(x, y, z, t) varies in both time and in some three-dimensional region with volume V. The verbal expression of conservation (7.1) remains valid. At any time t, the total amount of C in the volume can be computed by integrating c(x, y, z, t) over the volume total amount of C c(x, y, z, t) dV. (7.57) Now suppose that C is free to move about randomly, so that C moves in and out of the volume by passing through the volume's surface S. The flux J(x,...

## Diffusion in a Long Dendrite

All of the above examples examined steady behavior, after initial transients have decayed. However, reaction diffusion equations also contain information about the temporal evolution of the process to steady state. Consider calcium diffusing in a long dendrite. Suppose caged calcium is photore-leased from a small region around x 0. If we denote by c(x, t) the concentration of calcium along the length of the dendrite at each time t, then the model becomes dt Ddx ' < x < t > 0, (7.36) where...

## A41 Stability of Steady States

As we saw in the preceding section, a steady state may or may not be an attractor for nearby trajectories i.e., just because an initial condition is close to the steady state, it does not mean that after a time the trajectory will approach the steady state. However, when this is the case, the steady state is said to be stable and attractive or asymptotically stable. Three qualitatively different behaviors near steady state are illustrated by the solutions of the linear ODEs in Figure A.2. The...

## Two Component Oscillators Based on Autocatalysis

First, we consider two minimal requirements for oscillations in chemical reaction systems. 1. If the network is absurdly simple, with only one time-varying component, then oscillations are, generally speaking, impossible. (Proof For x(t) to be periodic, dx dt must take on both positive and negative values at some values of x, which is impossible if f (x p) is a continuous single-valued function of x.) Thus, to understand biochemical oscillations, we must start with two-component networks,...

## A52 Bifurcation at a Pair of Imaginary Eigenvalues

Limit cycles and periodic solutions are extremely important in physiology. Thus, one is often interested in whether or not they occur in a given system. Unlike fixed points that can be found exactly or graphically, it is much more difficult to determine whether or not there are limit cycles in a system. There is one method that is arguably the best and perhaps only systematic method of finding parameters where there may be periodic solutions in any system of differential equations. The...

## Calcium Oscillations in the Bullfrog Sympathetic Ganglion Neuron

Bullfrog sympathetic ganglion neurons (BFSG) are excitable cells with a full complement of voltage-dependent ion channels (see Yamada et al., 1998 ). However, the Ca2+ oscillations modeled here are driven by nonlinearity in the ER, with the plasma membrane playing only a passive role. As early as 1976, Kuba and Nishi Kuba and Nishi, 1976 observed rhythmic hyperpolarizations of the resting membrane potential when the neurons were exposed to caffeine. See Figure 5.2. The BFSG is a good model...

## Activation and Inactivation Gates

Channels can be thought to have gates that regulate the permeability of the pore to ions, as illustrated schematically in Figure 2.5. These gates can be controlled by membrane potential, producing voltage gated channels by chemical ligands, producing ligand gated channels or by a combination of factors. In a series of experiments, Alan Hodkgin, Andrew Huxley, and others established experimentally the voltage dependence of ion conductances in the electrically excitable membrane of the squid...

## Biochemical Kinetics and Feedback

In general, a biochemical reaction network is a schematic diagram of boxes and arrows. Each box is a chemical species. Solid arrows represent chemical reactions producing or consuming a molecule. Dashed arrows represent control by one species of a chemical reaction involving other species. Regulatory signals may be activatory (barbed end) or inhibitory (blunt end). To each reaction in a biochemical network is associated a rate law (with accompanying kinetic parameters). Hence, the network...

## Traveling Wave in the Fitzhugh Nagumo Equations

As we have seen in earlier chapters, chemical reaction schemes in cell biology can be quite complicated, involving many species. Furthermore, some species may be free to move, while others are not. Models of nerve axons, for example, include both diffusing species (transmembrane potential) and nondiffusing variables (the ion-gating variables, because ion channels are embedded in the membrane and do not move on the millisecond time scale of an action potential). Similarly, the Ca2+ wave induced...

## Exercises

Show that the assumed form of the indicator dye fluorescence (8.2) and the equilibrium relation (8.1) imply the backcalculation formula (8.3). 2. Assuming (8.2) and (8.4), derive the backcalculation formula for a ratiometric indicator, (8.5). 3. Beginning with (8.13), finish the derivation of the rapid-buffer approximation by substituting equilibrium values for Ca2+ T and CaB , A r Ca2 + + Ca2+ i B T _ V2 (D Ca2+ + D Ca2+ i B T + . dt (JCa i + Ca2+ i + k) _ V Tc i + Db Ca2+ i + k) + and taking...

## Suggestions for Further Reading

Random Walks in Biology, Howard Berg. This is a lovely introductory book on diffusion processes in biology Berg, 1993 . Mathematical Problems in the Biological Sciences, S. Rubinow. Chapter 5 gives a nice introduction to diffusion processes Rubinow, 1973 . Diffusional mobility of golgi proteins in membranes of living cells, N.B. Cole, C.L. Smith, N. Sciaky, M. Terasaki, and M. Edidin. This paper gives an example of how diffusion coefficients are measured in a specific biological context Cole...

## Positive Feedback Loop and the Routh Hurwitz Theorem

First, let us see whether Hopf bifurcations are possible in a system with a pure positivefeedback loop, i.e., a Jacobian of the form J+ at the steady state. Recall that the condition for a Hopf bifurcation is that the Jacobian matrix has a pair of purely imaginary eigenvalues A iw. The eigenvalues of J+ are roots of the characteristic equation G(A) A3 + (a + P + y)A2 + (aP + Py + Ya)A + - c A3 + AA2 + BA + C 0. The roots of this equation can be characterized by the Routh-Hurwitz theorem. Let...

## Negative Feedback Oscillations

The Jacobian J_ determines the stability of the steady state in a three-variable system with a pure negative feedback loop. In this case, G(A) A3 + (a + P + y)A2 + (aP + Py + Ya)A + aPY + c1c2 0, and the Routh-Hurwitz theorem implies that the steady state is unstable if and only if (a + P + 7)(aP + Py + 7a) < aPY + c . Furthermore, if equality holds, then G(A) (A + A)(A2 + B), so G(A) has conjugate roots on the imaginary axis at A aP + Py + 7a.

## Modeling Molecular Motions

The botanist Robert Brown first observed Brownian movement in 1827. While studying a droplet of water under a microscope he noticed tiny specks of plant pollen dancing around. Brown first guessed, and later proved, that these were not motile, although at the time he had no clue as to the mechanism of their motion. It was not until Einstein contemplated the phenomenon 75 years later that a quantitative explanation emerged. In order to develop an intuition about molecular dynamics we begin with...

## The Statistical Behavior of Polymer Growth

Since the intervals of time between events of monomers assembly and or disassembly are random, one can measure only the statistical behavior of polymer growth, such as the average velocity of the polymer's tip, (V) L(N) t, where L is the size of a monomer. Much useful information is buried in the statistical fluctuations about this mean velocity. One quantity that can be monitored as the polymer grows is the variance of the tip's displacement about the mean Var x(t) (x2) (x)2 L2((N2) (N)2). It...

## The Ensemble Density Approach Applied to the Stochastic Morris Lecar Model

The ensemble density approach described in Section 11.4.2 can be applied to the stochastic Morris Lecar model described above. The evolution equations for the conditional PDFs take the form d d Pa (n,v,t) - Jo (n,v,t) + a(v)(N - n + 1)PO(n - 1,v,t) (11.51) - a(v)(N - n)PO(n, v, t) + fl(v)(n + 1)Po(n + 1,v,t) - (v)nPa(n,v,t). This form is similar to (11.48) except that the transition probabilities are now voltage-dependent, and the probability fluxes JO(n, v, t) are given by the Morris Lecar...

## The Inositol Trisphosphate IP3 receptor

The IP 3 receptor is a Ca2+ channel located in the endoplasmic reticulum that is regulated both by IP3 and by Ca2+. Each receptor consists of three independent subunits, each of which must be in the open state for the channel to be open. Each subunit has one binding site for IP3 and two binding sites for calcium. Thus there are eight possible states for the subunit. Binding with IP3 potentiates the subunit. The two calcium binding sites activate and inactivate the subunit, and a subunit is in...

## Fluctuations in Macroscopic Currents

When the voltage clamp technique is applied to isolated membrane patches, openings and closings of single ion channels can be observed. Recall the single-channel recordings of T-type Ca2+ currents shown in the top panel of Figure 11.2. Importantly, the bottom panel of Figure 11.2 shows that when several hundred single-channel recordings are summed, the kinetics of rapid activation and slower inactivation of the T-type Ca2+ current are evident. In this summed trace, the relative size of the...

## Langevin Formulation for the Stochastic Morris Lecar Model

To consider the behavior of the Morris Lecar model under the influence of channel noise from a large number of K + channels, it is most convenient to use the Langevin formulation presented in Section 11.4.1. We do this by supplementing the rate equation for the average fraction of open K + channels, (11.50), with a rapidly varying forcing term where w is a random variable, ( ) 0, and the autocorrelation function of is given by h (t) (t')> Y(w, V)6 (t - t'). Following (11.36), y (w, V) is...

## The Resting Membrane Potential

The Nernst potential is the equilibrium potential for one permeant ion. In reality, no channel is perfectly selective for a given ion, and there are various channels selective for various ions in a given cell as well. The Goldman-Hodgkin-Katz (GHK) equation is related to the Nernst equation, but considers the case where there are multiple conductances. The GHK equation determines the resting membrane potential of a cell from a weighted sum of the various conductances V RT , Px K + out + PNa Na+...

## A1 Matrix and Vector Manipulation

Matrices can be multiplied, added, multiplied by scalar numbers, and differentiated according to the rules of linear algebra. Here we summarize these results for the 2 X 2 matrices and two-component column vectors . The matrix A multiplying the vector x acts as a linear operator that produces a new vector, z, according to the formula an a xi anxi + a x2 z Ax + . (A.1) a21 a22 j x2 j a21x1 + a22x2 j It can be verified using (A.1) that the identity matrix I 0 J leaves vectors unchanged, i.e., z x...

## Solving and Analyzing Differential Equations

Many students have worked with differential equations in their studies of the physical sciences or elementary mathematics, and they may have been introduced to solution techniques explicitly in an advanced calculus course. In general, the differential equations that arise for the rate of change in cellular properties will be complicated and difficult or impossible to solve exactly using analytical techniques. One example that Figure 1.8 A selection from the family of solutions to (1.6) for t 2....

## Htt 7 t

By methods of statistical physics beyond the scope of this book, the constant y( 0) can be shown to be inversely proportional to N and proportional to the sum of the rates of both the O C and C O transitions, that is, y ( o ) N- (n-36) An appropriate choice for is thus AB, where the AB are the increments of a Wiener process. Although we haven't fully justified this choice for , we can check that this random variable and (11.35) define the random variable fO in a manner consistent with the work...

## Ddix

J deprotonation at acidic reservoir + deprotonation at basic reservoir net protonation at both reservoirs Protonation rates are proportional to hyrogen ion concentration on either side of the membrane, and the net protonation at both reserviors is kp kp cacld + kpcbase. Note that the rates of protonation and deprotonation have the dimensions fcp (1 s) and kd (nm s), respectively. We assume that the deprotonation rates are the same at both reservoirs. First we nondimensionalize the model...

## Single Channel Gating and a Two State Model

The time course of voltage changes in a whole cell is the result of the average behavior of many individual channels. Our understanding of individual channel gating comes largely from experiments using the patch clamp technique (see Figure 11.1). For example, typical measurements from a so-called on-cell patch of T-type calcium currents in guinea pig cardiac ventricular myocytes are shown in the middle panel of Figure 11.2. The small current deviations in the negative direction indicate the...

## B15 Saving and Printing Plots

XPPAUT does not directly send a picture to your printer. Rather, it creates a PostScript file that you can send to your printer. If you do not have PostScript capabilities, then you probably will have to use the alternative method of getting hard copy outlined below. (Note that Microsoft Word supports the import of PostScript and Encapsulated PostScript, but can only print such pictures to a PostScript printer. You can download a rather large program for Windows called GhostView which enables...

## The Role of Mathematics

The scope of mathematical techniques employed to investigate problems in mathematical biology spans almost all of applied mathematics. The modeling of processes is discussed in detail here, but only the basics of the mathematics and the elementary tools for the analysis of these models are introduced. Rigorous mathematics plays at least three important roles in the computational modeling of cell biology. One role is in the development of the techniques and algorithms that make up the tools of...

## Excitability and Action Potentials

Another dynamical feature of the Morris Lecar model is excitability. A useful working definition of excitability is that a system is excitable when small perturbations return to the steady state, but larger (i.e., above a threshold) perturbations cause large transient deviations away from the steady state. An example of this is shown in Figure 2.11A, where the time course of the voltage with an applied current of 60 pA and initial conditions of w(0) wss 0.070 and either V(0) -22 or -17 mV are...

## Stability Analysis

The heart of answer lies in the fact that fixed points can be either stable or unstable and that the intersection of nullclines obtained when Iapp 150 pA is an unstable fixed point. Notice that the two trajectories that start nearest to the intersection of the nullclines in Figure 2.9B diverge away from the intersection and toward the limit cycle. The trajectory of the system is thus driven away from the intersection of the nullclines while at the same time being constrained to orbit around it....

## Molecular Mechanisms of Cell Cycle Control

Cell cycle events are controlled by a network of molecular signals, whose central components are cyclin-dependent protein kinases (Cdks). Cdks, when paired with suitable cyclin partners, phosphorylate many target proteins involved in cell cycle events (Figure 10.2A). For instance, by phosphorylating proteins bound to chromosomes at origins of replication (specific nucleotide sequences, where DNA replication can start), Cdks Figure 10.2 Cyclin-dependent kinase. (A) The role of a cyclin-dependent...

## The Goodwin Oscillator

The quintessential example of a biochemical oscillator based on negative feedback alone (Figure 9.3C) was invented by Brian Goodwin ( Goodwin, 1965, Goodwin, 1966 see also Griffith, 1968a ). The kinetic equations describing this mechanism are Vi Xi - k2 X2 , V2 X2 - As Xs . Here Xi , X2 , and X3 are concentrations of mRNA, protein, and end product, respectively v0, v1, and v2 determine the rates of transcription, translation, and catalysis k1, k2, and k3 are rate constants for degradation of...

## A42 Stability of a Nonlinear Steady State

What we have learned about stability of steady states for linear systems can be transferred partially to nonlinear ODEs. To be specific, let us consider a biological membrane with a gated ion channel. To do this we combine a model of ion gating with an expression that governs the membrane potential (see Chapter 1 and Chapter 2). For simplicity, we will consider only one conductance. If n represents the gating variable and V the voltage, then the two are coupled by the differential equations...

## Why Do Oscillations Occur

If the following three conditions on the Morris Lecar equations hold, then oscillations will occur the V nullcline has the inverted N shape like that in Figure 2.9B a single intersection of the V- and w-nullclines occurs between the maximum and Figure 2.10 (A) The K+, Ca2+, and total current (1Ca + 1K + ieak Iapp) when w 0.35. States 1 and 3 are stable steady states, and state 2 is unstable as indicated by the velocity vectors. (B) The phase plane for the Morris-Lecar model for Iapp 150 pA,...

## B22 Nullclines and Fixed Points

As discussed in earlier chapters, a powerful technique for the analysis of planar differential equations and related to the direction fields is the use of nullclines. Nullclines are curves in the plane along which the rate of change of one or the other variable is zero. The x-nullcline is the curve where dx dt 0, that is, f (x, y) 0. Similarly, the Figure B.9 Direction fields and some trajectories for the FitzHugh-Nagumo equations. y-nullcline is the curve where g(x, y) 0. The usefulness of...

## A4 Phase Plane Analysis

Obtaining a solution to first-order ODEs means that you have expressed all of the dependent variables as functions of time. In the case of the 2 X 2 linear equations in Section A.3, this means that we have the time series for x, and x2. A great deal can be learned about these solutions by plotting the dependent variables as a function of time as done in Figure A.1. However, there are other ways of plotting solutions that give additional insight. For example, one can plot x, versus time, or some...

## The Voltage Clamp

In order to measure the voltage across a cell membrane or the current flowing through a membrane, microelectrodes are inserted into cells. These electrodes can be used both to measure current and voltage and to apply external current. In order to measure the Figure 2.7 Simulation of voltage clamp experiment using (2.28) and (2.29). (A) Current records resulting from 40 ms depolarizations from the holding potential of -60 mV to the indicated test potentials. (B) the maximum (steady state)...

## B17 Looking at the Numbers The Data Viewer

In addition to the graphs that XPPAUT produces, it also gives you access to the actual numerical values from the simulation. The Data Viewer shown in Figure B.8 has many buttons, some of which we will use later in the book. The main use of this is to look at the actual numbers from a simulation. The independent variable occupies the leftmost column, and the dependent variables fill in the remaining windows. Click on the top of the Data Viewer to make it the active window. The arrow keys and the...

## Traveling Wave Solutions

An interesting and important problem is to determine the behavior of the bistable equation when a portion of the region is initially above the threshold a and the remainder is initially at zero. To get some idea of what to expect it is useful to perform a numerical simulation. For this numerical simulation we use the method of lines to solve the differential equations ddt0 ax (ci(t) - c (t))+f (c ), (7.77) dCn D (c +i (t) - 2c (t) + c _i (t)) + f (c ), n 1,2, ,N - 1, (7.78) (cn-i(t) - Cn(t))+ f...

## Glucose Dependent Insulin Secretion

Insulin is secreted from -cells in the pancreas in an oscillatory fashion. Glucose must be metabolized by the -cell to stimulate insulin secretion, and the insulin, which is Figure 4.4 Flow system for experimental study of insulin secretion. Vbed is the volume of the reaction bed, f is the volume flow rate, and Go is the inflow concentration of glucose. Redrawn from Maki and Keizer, 1995 . prepackaged in secretory vesicles, is secreted from the -cell into the capillary system by exocytosis....

## John J Tyson

Biochemical and biophysical rhythms are ubiquitous characteristics of living organisms, from rapid oscillations of membrane potential in nerve cells to slow cycles of ovulation in mammals. One of the first biochemical oscillations to be discovered was the periodic conversion of sugar to alcohol glycolysis in anaerobic yeast cultures Chance et al., 1973 . The oscillation can be observed as periodic changes in fluorescence from an essential intermediate, NADH see Figure 9.1. In the laboratories...

## Phase Plane Analysis

The mechanistic features underlying continuous spiking and action potentials can be understood easily using phase plane analysis. Because of the ease of representation in two dimensions, two variable models such as the Morris Lecar model are particularly amenable to this technique. Phase plane analysis is a powerful way to determine how the behavior of a system will change with changes in the various parameters in a system. Several types of plots are utilized as part of what is generically...

## Equations for Membrane Electrical Behavior

The Nernst potential represents the chemical equilibrium, and the GHK equation represents the steady state. Given the several conductances and their reversal potentials, we calculate what the membrane potential will be after the system has stabilized. These equations tell us nothing about how the system evolves to the steady state. Because we are interested in the time course of the membrane voltage, we have to study the dynamics of the various currents that flow in and out of the cell. We can...

## Model of the Fertilization Calcium Wave

When mature Xenopus laevis oocytes eggs are loaded with an indicator dye e.g., Ca2 -green dextran and stimulated by the fusion of sperm, a propagating wave of intracellular Ca2 release can be observed by backcalculating the free Ca2 concentration Ca2 i x, t from a time-dependent fluorescence signal F x,t according to 8.3 or 8.5 . This fertilization Ca2 wave is an important step in early development. It triggers Figure 8.2 Schematic diagram of the fertilization Ca2 wave model. Ca2 enters the...

## Numerical Integration of Differential Equations

Even if 1.14 were more complicated and could not be solved exactly, a numerical approximation could still be calculated. The simplest and perhaps the oldest method of numerical solution goes back to the mathematician Euler and is easy to understand. The method is called the forward Euler method and it is a prototype for all other Figure 1.9 The Euler method of numerical integration relies upon a series of short linear approximations using the derivative at the old time point. The solution to...

## Flux Balance Equations with Rapid Buffering

We begin with a general description applicable to any cell with a cytosolic compartment subscripted by i, originally for intracellular and an ER subscripted by ER and one species of Ca2 buffer in each compartment. That gives a total of four species, with bound and free Ca2 in each of the two compartments. We define N total number of Ca2 ions, bound and free, in the cytoplasm and NER total number of Ca2 ions in the ER. Then Ca2 r N, Ca2 EoR ER, 5.1 where the V's representing the volumes of the...