Barnes D.J., Chu D. Introduction to Modeling for Biosciences

Подождите немного. Документ загружается.

5.6 Maxima: Case Studies 213

parameters. We decide to keep the parameters K and h ﬁxed, and only vary b.

For this we deﬁne the function genres that takes as inputs our choice for these

parameters, and then generates results for various values of b.Todothisweusea

for ... step ... thru ... do construct, that allows us to deﬁne loops

in Maxima. In the genres function, i is an index variable that simply counts where

we are in the loop at a particular time. During each loop it will be incremented by

a value speciﬁed after the step keyword; in this case we chose it to be 1/100. It

then executes a command (or set of commands) as speciﬁed after the do keyword.

Notice that the commands in the loop’s body must be enclosed within parentheses.

The loop stops when the index variable i reaches the number speciﬁed after the

thru keyword; in this case 1/2.

The task of this loop is to generate results for various values of b. An inspec-

tion of the function deﬁnition below shows that the value of i is used as a value

for b. The inner loop in genres (with index variable j) puts the results from

ssv into a list called finL. Each member of the list is itself a list of the shape

[bvalu,ssvalu], where bvalu is the value of b (or i) and ssvalu is a steady

state value for this particular parameter.

(%i5) genres(pK,ph):= block( finL: [[0,0]],

for i:1/10 step 1/100 thru 1/2 do (

tmpres: ssv(pK,ph,i),

for j:1 step 1 thru length(tmpres) do(

finL:append(finL,[[i,rhs(tmpres[j])]]))),

finL: delete([0,0],finL), finL)$

(%i6) finalV: genres(2,3)$

(%i7) plot2d([discrete,finalV],[style,[points]]);

In %i6 we call genres, which returns the list of steady states in convenient

shape to be plotted by plot2d as discrete points. We do not actually show the

resulting plot (but see Fig. 4.17 on p. 178); the reader who has actually typed in our

code will notice that the command on %i8 would also plot some negative steady

states, which are, of course, undesired. The simplest solution is to restrict the range

of the plot to positive values. Finally, we determine the stability properties at b

1/3. First we generate the Jacobian, then we substitute the steady state values for the

particular set of parameters into the Jacobian to be able to determine the stability.

We know from the bifurcation diagram (Fig. 4.17 on p. 178) that at b =1/3 only

one steady state exists.

(%i8) load(linearalgebra)$

(%i9) J: jacobian([eq1,eq2],[x,y])

(%o9)

⎛

⎝

−b −

h−1

(

)

−

h−1

(

)

−b

⎞

⎠

(%i10) Jv: ev(J,K=2,h=2,b=1/3)

(%o10)

⎛

⎝

−

(

)

−

(

)

−

⎞

⎠

214 5 Mathematical Tools

(%i11) eigenvalues(ev(Jv,x=rhs(ssv(2,2,1/3)[1]),

y=rhs(ssv(2,2,1/3)[1]) ))$

(%i12) float(%)

(%o12)

[[

−.6172337335394495, −.04943293312721716

]

[

1.0, 1.0

]]

Line %o12 shows the two eigenvalues and their multiplicities. Since both eigen-

values are negative, we conclude that the steady state point is stable.

5.7 Summary

Maxima is a very powerful computer algebra system. Compared to many other sys-

tems, it is intuitive to learn, easy to use, fast and portable. Best of all, it is freely

available.

In this chapter we have provided an introduction to Maxima. This is enough to

solve the problems discussed earlier in this book and to get started on an independent

exploration of this system. The sheer breadth of Maxima means that only the most

important commands and concepts have been discussed here and much has been

left out. While using the system, the reader will nearly certainly encounter problems

and be “stuck” for some time. In our experience, Maxima is very user friendly and

contains few of the idiosyncrasies that plague other systems. Nonetheless, being

stuck is part of the research process, as frustrating as it seems!

There are a number of web-based resources to help the reader. One of the best

points of departure is the Maxima home page

which provides a number of links to

tutorials and manuals of varying lengths, depths and qualities. Particularly valuable

is the Maxima reference manual that lists all the Maxima commands and packages

with a brief description.

An aspect that has been nearly completely left out here is the ability to program

Maxima. Most users will probably not wish to write extensions to Maxima in Lisp.

On the other hand, most will, sooner or later, write more complex procedures using

the Maxima command language. In a rudimentary way we have already given an ex-

ample above, but it is not difﬁcult to write more advanced and powerful procedures

and programs in Maxima.

http://maxima.sourceforge.net.

Chapter 6

Other Stochastic Methods and Prism

If we model a system using differential equations then we make a number of implicit

assumptions about the nature of the system. These assumptions will sometimes be

correct, sometimes incorrect but reasonable, and sometimes incorrect and unreason-

able. At the heart of the method of differential equations is the idea of inﬁnitesimal

changes. Remember that a statement of the form,

˙x =f(x)

means that the state variable x changes during every inﬁnitely small window of time

dt like f(x)dt. Assume now that our differential equation describes the change in

concentration of a molecular species due to chemical reactions. The assumption un-

derlying any such differential equation is that the concentration is a continuously

varying quantity; so it could take a value of 0.3012 or 0.4919 or anything in be-

tween, above or below. Fundamentally, of course, chemical systems are made up of

individual particles. Concentrations are just the number of particles per volume and

are, therefore, only approximately continuous. This means that there is discrepancy

between the assumptions made in deterministic models and reality. This discrepancy

by necessity leads to an error of the model whenever nature happens to be discrete.

It clearly makes no sense to speak of a change of the number of molecules by

0.001 or, indeed, any other non-integer change of particle numbers. Yet it is pre-

cisely such a fractional change that a differential equation model may describe. In-

deed, as the system approaches the steady state, changes in the particle numbers will

become successively smaller, and eventually require changes of concentrations that

correspond to fractions of molecules.

To see the origin of the modeling error, consider a system where molecules of

type A are converted into molecules of type B. Assume that a differential equation

model predicts that, at time 10, 4.7361 molecules of A will be converted into B.

Being discrete, the real system can only change by either exactly 4 or 5 molecules.

Which one is the right solution? The value from the differential equation gives us

a hint that the conversion of 5 molecules is more likely, but we cannot be sure.

In either case, even if up to t = 10 the differential equation model has accurately

described our system (which it most likely has not) then at least at time t =10 +t

D.J. Barnes, D. Chu, Introduction to Modeling for Biosciences,

DOI 10.1007/978-1-84996-326-8_6, © Springer-Verlag London Limited 2010

215

216 6 Other Stochastic Methods and Prism

the prediction will necessarily be at odds with reality. Our assumed natural system

can, by its very nature, not follow the differential equation. Depending on whether

4 or 5 molecules have been converted, the differential equation will predict a higher

or lower concentration of A molecules than exhibited in the real system.

When we choose to model the change of discrete entities (such as molecules) us-

ing a formalism that is inherently continuous in nature (differential calculus), then it

is clear that there will be modeling errors. This tells us that the differential equation

models are always inaccurate models when they describe discrete systems. The rel-

ative importance of that inaccuracy will depend on the size of the system. A relevant

measure could be the number of particles in the system. The higher this number, the

better the deterministic models. For very small systems, deterministic approaches

are often too inaccurate to be useful.

This problem is further exacerbated by truly stochastic effects in biological sys-

tems. More often than not, real systems do not follow a deterministic law but their

behavior has elements of randomness. On average, these systems may often follow

a deterministic law and are accurately described by differential equations. Funda-

mentally though, their apparently deterministic behavior is the emergent result of

a myriad of random effects. In very large systems, the stochastic nature of these

systems may not be noticeable. For smaller systems, however, these ﬂuctuations

become increasingly dominant. The actual behavior of these small systems will sig-

niﬁcantly deviate from their mean behavior. This deviation of the actual behavior of

a system from the mean behavior is often referred to as noise, and is the result of

stochastic or statistical ﬂuctuations. In the literature, these words are largely used

interchangeably to describe one and the same phenomenon.

This means that small and discrete systems with stochastic ﬂuctuations are typi-

cally not very well described by differential equation models. There are, of course,

huge advantages to be gained from modeling systems using differential equations.

Theoretically, this method is very well understood and, often, statements of con-

siderable generality can be deduced from differential equation models. Moreover,

computationally they are relatively cheap to analyze numerically, particularly when

compared to explicit simulation methods. Altogether, there are signiﬁcant argu-

ments in favor of using differential equations to model systems. Hence, even if a

system has signiﬁcant stochastic effects and only approximately follows a determin-

istic path, one may still prefer differential equations to explicit stochastic models. In

modeling, the question is never whether or not a model is wrong, but always, how

useful it is.

If our system consists of millions of molecules that interact with one another,

then it will still be a system consisting of discrete entities. At the same time it will,

for all practical purposes, behave very much like a continuous system. In such cases,

the differential equation approaches will remain very accurate. The problem starts

when the number of particles is not very high, and when the discrete nature of the

particles starts to manifest itself in the phenomenology of the system. In biology,

such cases are common. Transcription factors, for example, are often present in

copy numbers of thousands or maybe only hundreds of molecules. In these cases,

the deviation of the actual behavior from the deterministic behavior predicted by the

differential equation can be quite substantial.

6.1 The Master Equation 217

An example of a system where random ﬂuctuations are often taken into account

is gene expression. There are many layers of stochasticity in the system. For one,

the gene may be controlled by a stochastic signal (for example, a ﬂuctuating tran-

scription factor concentration). Secondly, binding of any regulatory molecules to

their sites is stochastic. Thirdly, even under constant external conditions, the ac-

tual expression of the protein is a stochastic process, in the sense that the number

of proteins produced will ﬂuctuate from one time interval to the next. Finally, the

breakdown or loss of proteins from the cell is also a stochastic process in this sense.

Taking into account that the total molecule number may be very low (hundreds of

molecules rather than millions) the ﬂuctuations of the protein concentration may

dominate the system and cannot be ignored any longer.

Protein expression is but one example of an area in biology where stochastic ef-

fects are often important. There are many more. Noise and stochastic ﬂuctuations

occur at all levels of organization in biology and their importance is increasingly

recognized. Consequently, over recent years signiﬁcant research effort has gone into

both developing new methods and tools to model noise, and to applying those meth-

ods to biological problems. Much progress has been made in understanding and

describing stochastic systems in general. As a result, there are a number of mathe-

matical techniques, as well as computational tools and algorithms, available to deal

with stochastic systems.

The focus of this chapter is to give an introduction into mathematical techniques

to deal with stochasticity. As (nearly) always, mathematical approaches are prefer-

able to simulations. Unfortunately, when it comes to stochastic systems the math-

ematics involved is rather difﬁcult and even relatively basic models require a level

of mathematical skill that is normally the reserve of the expert. Moreover, as with

most mathematical modeling methods, realistically sized problems quickly test the

limits of tractability and will challenge even the skills of the specialist. Given the

importance of these approaches it is nonetheless useful to understand the basics of

stochastic modeling techniques, as it will provide some essential background mate-

rial to understand the primary literature in the ﬁeld.

Within the scope of this book it is only possible to give a taster of what can be

done and what has been done mathematically. The reader who wishes to deepen her

mastery of the subject will ﬁnd Gardiner’s “Handbook of Stochastic Methods” [18]

a useful starting point.

6.1 The Master Equation

The basis of many mathematical models of stochastic systems is the so-called mas-

ter equation. The name “master equation” is a bit misleading. In essence, it is a

differential equation, or set of differential equations, that formulates how the proba-

bilities associated with a particular system change over time. One way to think about

it is to assume that a system can be in a number of states, s = 1, 2, 3,...,n, each

with a certain probability P

. Concretely, the P

could represent the probability that

there are s proteins in the cell, or that s molecules are bound to the DNA, and so on.

218 6 Other Stochastic Methods and Prism

The numbers here are just labels and should not be understood to refer directly to

some aspect of the model. We have here assumed that there is only a ﬁnite number

of such states, but one could equally imagine that the index runs from −∞ to +∞.

The master equation is then a set of differential equations, each of which de-

scribes how the probability of being in a particular state s changes over time.

,...,P

) (6.1)

As is usual, we have suppressed the time dependence of the probabilities, and we

have also assumed that the states are labeled from 1 to n. In principle, the solutions

to the master equation give us the probabilities of being in the various states as a

function of time, which is the information one normally wishes to have in stochastic

modeling applications. The problem is that, for the moment we do not know the

function f

; we will describe below how this can be found. Another problem is

more serious: Even when f

is known, in most cases it is not possible to solve

the master equation exactly. In some cases good approximations can be found, but

this is nearly always difﬁcult. In the vast majority of cases one needs to resort to

simulations to solve the master equation. Yet, despite being rarely of great utility, it

is still worth being familiar with the concept of the master equation because it plays

a very important role in stochastic modeling in biology.

One way to think of the master equation is to regard it as a gain-loss equation of

probability similar to the normal differential equations encountered in Chap. 4.An

example to illustrate this idea is the random walk in 1 dimension. The basic setting

is that there is an inﬁnite number of squares lined up such that each square has

exactly two neighbors. The “random walker” can, at any one time, occupy exactly

one of these squares. From its current position the walker can move to one of the two

adjacent sites with a rate of k, or choose to remain where it is. (We assume transition

rates and hence continuous time. The same problem could also be formulated in

discrete time.) The relevant states for this problem are the occupied squares, that is

we can assume the system to be in state s at time t if the random walker is on the

square labeled s at time t . Since we assumed an inﬁnite number of squares, for this

particular problem it is then also most suitable to let the state labels run from −∞

to +∞.

To formulate the master equation for this problem it is helpful to ﬁrst remember

what exactly we want to know: In the absence of any information, it is impossible

to make any statements about the whereabouts of the random walker, except that it

could be anywhere. The situation changes when we start to take into account that, at

time t =0 the random walker was on square 0. Given this information/assumption,

what we are interested in is the probability that the random walker is at square 45,

say, at some point in time, say 10 minutes after starting. Mathematically, this is

often formulated using the notation from conditional probabilities. The expression

P(s,t|0, 0) denotes the probability that the random walker is on square s at time t

given that it was at square 0 at time 0.

We cannot write down the probability of this outright. However, what we can

do is to write down how this probability changes. If we assume that at time t the

random walker is actually at square s, then we know that, with a rate of k, it will

6.1 The Master Equation 219

move to the left or to the right. Altogether, it thus moves away with a rate of 2k.

If we express this in terms of probabilities, we can also say that the probability

of being on s reduces with a rate of 2k. Of course, the random walker does not

just disappear, but moves on to one of the adjacent sites. One man’s loss is another

man’s gain, so to speak. Our focal square s is the neighbor of the two other sites

and can accordingly proﬁt from their loss. From each of the adjacent sites it will

experience an inﬂux with a rate of k. Altogether we can thus formulate the change

of probabilities at each site in terms of the rates as a differential equation.

∂P(s,t|0, 0)

∂t

=kP(s −1,t|0, 0) +kP (s +1,t|0, 0) −2kP (s, t|0, 0) (6.2)

This is the master-equation for the random walk.

Often one sees a more compact way of writing the master equation using the jump

operator, E, which transforms expressions in terms of s into expressions in terms

of s +1; the inverse of this operator, E

−1

that transforms s into s − 1. The jump

operator does not correspond to any entity in real life, but is merely a convenient

abbreviation tool. We will not use it any further within this book but merely bring

it to the reader’s attention as she will almost certainly encounter it in the literature.

Using this operator we obtain a different form of the master equation (6.2).

∂P(s,t|0, 0)

∂t

=k(E

−1

+E −2)P (s, t|0, 0) (6.3)

6.1 Show how to derive (6.3) from (6.2).

Let us consider another simple example. Assume that the probability of an event

happening is proportional to some waiting time. For example, if we “sit” on a

molecule in a well stirred vessel, then the probability P

t

of colliding with an-

other molecule within the next time unit t is (more or less) proportional to t.

Using m as a proportionality constant, we can express this formally as

t

=mt (6.4)

Assume now that we are interested in the number of collisions since time t = 0.

Let us denote by P(k,t) the probability that at time t there have been exactly k

collisions between the molecule we are sitting on and others. We do not know this

probability but, given our assumption in (6.4), we can formulate how it changes

P(k,t+t) = P(k,t)(1 −P

t

) +P(k−1,t)P

t

= P(k,t)(1 −mt) +P(k−1,t)mt (6.5)

This equation can be interpreted very easily. It consists of the sum of: the probability

that there have already been exactly k collisions by time t (given by P(k,t)) and

that another does not follow within t (given by (1 −P

t

)); and the probability

that there have been k −1 by time t (given by P(k−1,t)) and that a single colli-

sion follows within the next t. One could imagine that within t there could be

220 6 Other Stochastic Methods and Prism

more than a single collision, which would complicate the calculation considerably.

A relatively easy way out of this problem is to assume that the time interval, t ,

is so short that more than a single collision is unlikely enough to be ignored as a

possibility. This is precisely what we do. We can now further transform (6.5).

P(k,t+t) −P(k,t)

t

(

P(k−1,t)−P(k,t)

)

(6.6)

Taking the limit t → 0 we obtain a differential equation for the probability of

observing k collisions within a time period of t.

∂P(k,t)

∂t

(

P(k−1,t)−P(k,t)

)

(6.7)

This is a master equation. At a ﬁrst glance, it seems very similar to the simplest

differential equations we have encountered in Chap. 4 and we would expect it to be

easily solvable using elementary methods. However, a closer examination reveals

a complication. The ﬁrst argument to P on the right-hand side of the equation is

not the same as that on the left-hand side. The symbols P on the right hand side

are therefore different functions of t. Since the left-hand side of the differential

equation and the right-hand side do not contain the same functions, our approaches

from Chap. 4 are not applicable, at least not immediately.

To solve this master-equation we must resort to a trick, using the so-called prob-

ability generating function.

G(x, t)

∞



k=0

P(k,t) (6.8)

This function is a transformation of the probability functions P . The reader will

notice the new variable x that the transformation introduces, while dispensing with

the original variable k (over which the sum is taken). The easiest way to think of

x is not to think about it at all, and just take it as a formal variable without any

intrinsic meaning! As it will turn out, the master equation (6.7) will be very easy

to solve for this transformed function. Once a solution is obtained, we can then

transform back to obtain the desired solution. This approach works because the

generating function sums over the offending variable k of the (as yet unknown)

probability function of observing exactly k collisions. At the same time, it leaves

the time variable, t untouched. Since x is not dependent on time, there is a clear

relationship between the time derivative of G and the probability function.

∂G(x,t)

∂t

∞



k=0

∂P(k,t)

∂t

This expression is simply obtained by differentiating (6.8) with respect to time.

We can now reformulate the differential equation (6.7) in terms of the probability

6.1 The Master Equation 221

generating function. To do this, we ﬁrst multiply (6.7)byx

, and then take the sum

over all k.

∞



k=0

∂P(k,t)

∂t

∞



k=0



P(k−1,t)−x

P(k,t)



= m



∞



k=0

P(k−1,t)−

∞



k=0

P(k,t)



= m



∞



k=0

P(k,t)−G(x, t)



∂G(x,t)

∂t

= m(x −1)G(x, t)

(6.9)

We can now solve (6.9) using the methods of Chap. 4.

G(x, t) =G(x, 0) exp

(

m(x −1)t

)

(6.10)

The ﬁrst question to address is to determine the value of G(x, 0) =



P(k,0).

Given the meaning of P we must assume that at time t = 0 there have been

no collisions at all yet; we have, so to speak, reset the counter at t = 0. Hence,

G(x, 0) =P(1, 0) =1.

Now, we have a solution for the probability generating function, but not yet for

the probability density itself. To obtain that, we will use the Taylor series for e



/n! to expand the right hand side of (6.9) around 0.

G(x, t) = exp

(

mt (x −1)

)



P(k,t)

exp(mt)

(mt) x

exp(mt)

(mt)

2! exp(mt)

(mt)

3! exp(mt)

(mt)

4! exp(mt)

+···

= exp(−mt)



(mt)

(6.11)

The solution can now be obtained by comparing the coefﬁcients of the expansion of

G(x, t) with the coefﬁcients of the sum



P(k,t). We then obtain the probability

function for every k.

P(k,t)=

(mt)

exp(−mt)

(6.12)

The reader may recognize this expression as the Poisson distribution.

6.2 Fill in the missing steps to get from (6.11)to(6.12).

A similar method can be used to solve the master equation for the random walk

(6.2) above. Remember that we denoted by P(n,t|0, 0) the probability of the ran-

222 6 Other Stochastic Methods and Prism

dom walker being at n at time t given that he started at 0 at time t =0.

∂P(n,t|0, 0)

∂t

=kP(n −1,t|0,0) +kP(n +1,t|0, 0) −2kP(n,t|0, 0) (6.2)

This equation can be solved using the characteristic function, which is similar to the

probability generating function.

G(x, t) =exp(ixn)=

∞



n=−∞

exp

(

ixn

)

P(n,t) (6.13)

where i is the complex unit, that is the square root of −1. Note that here we have

omitted the initial conditions on the probability to simplify the notation. As before,

we can transform the differential equation (6.2) into a differential equation of the

characteristic function by multiplying both sides with exp(ixn) and summing over

all n. Following a similar procedure as in the derivation of (6.9), we obtain a new

equation.

∂G(x,t)

∂t

=k(exp(ix) +exp(−ix) −2)G(x, t) (6.14)

6.3 Fill in the details of how to reach (6.14).

Again, this can be solved for G(x, t) using the methods from Chap. 4.

G(x, t) =G(x, 0) exp



(exp(ix) +exp(−ix) −2)kt



(6.15)

As in the ﬁrst example, we can now determine that G(x, 0) = 1 by choosing our

initial conditions properly. Unfortunately, this leaves us with an unpleasant expo-

nential that has an exponential as its argument. What makes the situation worse is

that this argument is not even real, but complex. However, things are not as bad as

they seem. To solve this we can use the fact that exp(±ix) =cos(x) ±i sin(x).

G(x, t) =exp[(2cos(x) −2)kt] (6.16)

To further process this, we expand the cosine around x =0 using a Taylor series.

cos(x) =1 −

+··· (6.17)

We next make an approximation to keep the problem simple. Instead of keeping the

entire Taylor expansion, we will only keep terms of order x

and lower. To truncate

Taylor expansions in this way is a very common technique in physics to simplify

expressions. It usually means making a mistake, in that the truncated expansion is

only an approximation to the original function. One hopes that this mistake is very

small and does not matter in practice. In any case, a small inaccuracy is mightily

outweighed by the massive simpliﬁcation gained.