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

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

6.5 Analyzing Markov Chains: Using PRISM 263

Fig. 6.4 The property speciﬁcation interface

6.5.5 The PRISM GUI

PRISM can be conveniently used from the command line but also has a graphical

user interface, including an editor and a facility to check properties and plot results

from several model checking queries. The graphical user interface can be called with

the xprism command. The interface is reasonably intuitive to use and the reader

should be able to quickly ﬁnd her way around without any instruction, so we omit a

detailed description here. Nevertheless, we include a few tips and tricks that may be

helpful.

The PRISM GUI consists of several tabs: the model tab containing the editor, the

properties tab with the model checking interface where properties can be speciﬁed,

analyzed and results plotted, a simulation tab where numerical sample paths can be

created, and a log tab that reports details about the Markov chain (see bottom left in

Fig. 6.4).

The advantage of the graphical user interface is that it allows the user to automate

querying of PRISM while varying parameter values. The results can then also be

plotted as a function of the variable parameter. PRISM more or less automatically

varies parameters that are left unspeciﬁed in the model deﬁnition ﬁle. The usage of

this feature is easy but not entirely intuitive. To illustrate the general principle of

interacting with the PRISM GUI we will walk the reader through one example. The

We are referring here to PRISM version 3.3.

264 6 Other Stochastic Methods and Prism

objective of the example is to plot the probability of b =1 as a function of time for

the model prModel3.pm.

The property checking interface of PRISM consists of four different panes, for

entering properties, entering constants, deﬁning experiments, and the plot pane that

contains all the plots (shown in Fig. 6.4).

1. In the editor tab, load the model ﬁle prModel3.pm into PRISM and switch to

the properties tab.

2. Move the mouse to the properties pane, right click, choose Add and type in the

property speciﬁcation, i.e., P=?[true U[T,T](b=1)]. Note that we leave T

unspeciﬁed here. PRISM will complain about this but the reader may rest assured

that ignoring this complaint will do no harm.

3. Next, right click on the contents pane and choose Add constant. In the name

ﬁeld type T, and in the drop down list under type choose double. This means

that we specify our constant T to be continuously varying. Leave the value un-

speciﬁed.

4. Go back to the properties pane and click onto the line containing the previously

entered property until the whole line becomes light blue. The properties pane can

take more than one property speciﬁcation and this step serves to specify which

property should be checked. The next step will not work unless this step has been

completed successfully and the line containing the query is light-blue.

5. In the experiments pane, right-click again and choose New experiment.This

brings up a window where can be speciﬁed the range over which T should be var-

ied and the step size. This is self-explanatory. The graph in Fig. 6.4 was produced

by varying T from 0 to 10 using a step size of 0.1.

6. Following this, a window will come up asking whether to plot the results into a

new or existing window. After conﬁrming this step the graph of all probabilities

will be plotted.

Generally, in order to produce a graph it is necessary to leave a constant unspeciﬁed.

This could be a part of the property speciﬁcation itself, or constants in the model

deﬁnition can also be left unspeciﬁed. Through experimentation, the reader will

quickly ﬁnd out how to do this. Right-clicking at the appropriate places also allows

the user to see the raw data used to produce the graph. The graph itself can also be

exported into various formats, including Encapsulated Postscript (EPS) and Matlab

format.

6.6 Examples

In the remainder of this chapter we provide two case studies of the use of PRISM to

analyses stochastic systems. The ﬁrst case study will be a detailed model of the ﬁm

switch from Chap. 2. The second case study will illustrate how the simulation func-

tion can be used to produce stochastic simulations when given a set of differential

equations.

6.6 Examples 265

6.6.1 Fim Switching

The ﬁm switch in E.coli is regulated by the invertible genetic element ﬁmS (see the

discussion in Sect. 2.6.2.2). There are two binding sites to each side of this element.

Altogether, there are two types of proteins that can bind to these competitively,

namely FimB and FimE. The element ﬁmS is inverted when all four binding sites are

occupied by the same type of molecule. “Inverting” literally means that the element

is cut out of the genome and re-inserted in the opposite orientation (Fig. 6.5). This

particular switch regulates the expression of so-called ﬁmbriae in E.coli; these are

small hair-like protrusions from the surface of the cell. Their primary role is to aid

attachment to other cells. For E.coli cells that live within the gut of other organisms,

this means that they are able to hold on to tissue from the host.

Fimbriae are virulence factors. In the context of commensal cells this means that

they need to be regulated to avoid the host from triggering its defenses which would

lead to an inhospitable environment for the parasites. Interestingly, the way nature

chose to regulate ﬁmbriation is not at a global level, for example via cell-cell com-

munication. Instead, ﬁmbriae are regulated as a stochastic function at the level of the

individual cell. Each bacterium switches randomly between the ﬁmbriate and aﬁm-

briate states. The levels of ﬁmbriation in the population are controlled indirectly

by the environmental conditions through the switching bias from aﬁmbriate to ﬁm-

briate. As a result, the population as a whole is a mix of ﬁmbriate and aﬁmbriate

cells.

Mechanistically, the switching rates to and from the ﬁmbriate state emerge from

the competition of two proteins, or recombinases, for the four binding sites near

ﬁmS. The recombinase FimB can switch ﬁmS in both directions, however, it does so

with a very low rate. The other recombinase, FimE only switches ﬁmS from the on

state (where ﬁmbriae are expressed) to the off state (where they are not). While it

switches only in one direction, it does so with a much higher efﬁciency than FimB.

FimE is only expressed when the ﬁmS switch is in the on orientation.

We will now present a simpliﬁed model of this system to demonstrate the prac-

tical use of PRISM. To keep the model transparent we will make a number of key

simpliﬁcations:

• The concentrations of FimB and FimE are constant over time.

• Binding to each of the binding sites takes place with a rate that is proportional to

the concentrations of the FimB and FimE proteins respectively.

• On each side of ﬁmS, if one of the two binding sites is occupied by, say FimB,

then FimE cannot bind to the other site. There are no such constraints between

the double sites on each side of ﬁmS.

In the model shown in Code 6.3, we have split the Markov chain over three modules:

one to keep track of the switching, and one each to keep track of the binding at each

site of ﬁmS. In the PRISM speciﬁcation code below, we will only explicitly write

one of the binding modules (the one for the “left” binding site). The right module

is exactly the same as the left module. We can therefore simply copy the existing

266 6 Other Stochastic Methods and Prism

Fig. 6.5 A schematic illustration of the ﬁm switch. FimB and FimE compete for binding sites. If

all four binding sites are occupied by one species, then the inversion of ﬁmS can happen. In the

schematic ﬁgure on the left, all binding sites are occupied and the switch is ready to invert. The

right-hand side shows the inverted switch. The binding sites are not fully occupied so no switch

can take place at this stage

module and rename the state variables. PRISM provides a function for this, and an

example of the syntax is given in Code 6.3.

We have omitted the speciﬁcation of the constant rate factors here (to save space).

The reader who wishes to test this model is encouraged to obtain a copy from the

book’s web site and experiment with her own values.

Biological modeling in PRISM can be very intuitive. In this case, we represent

each of the two binding sites in their own modules. The state variable bl counts how

many molecules of FimB are bound to the left binding site; el is the analogous state

variable for FimE on the left binding site. In the copied modules rightBS we have

the corresponding state variables br and er.

For each of the binding sites there are only three transitions, namely from no

occupation to one molecule being bound, and from one molecule bound to two

molecules bound. We also allow for unbinding of molecules. Note that the state

transitions disallow a mixed bound state, where on the left side, say, one FimB and

one FimE are bound.

The inversion of ﬁmS is represented in the switching module. The transitions

switch1 and switch2 formulate the switching on of the system. The former is

the more important, as it controls switching when all four binding sites are occu-

pied by FimB; this is represented by the conditions bl = br = 2. The latter repre-

sents switching ﬁmS from off to on when all four sites are occupied by FimE. We

know empirically that this is a very rare event and therefore one should set the rate-

constant eeswitchon to a very low value. The same module also represents the

switching off of ﬁmS. Note that the conditions for switching are such that switching

only occurs when all binding sites are occupied by the same type of protein.

Using this model (together with the parameters that can be found in the online

version of the model ﬁle) we asked PRISM about the steady state probability of the

switch being on, i.e., in state 1. Using the PRISM GUI we generated a graph (see

Fig. 6.6) that shows these probabilities for various concentrations of FimB.

6.6 Examples 267

module leftBS

ctmc

const double btobb = 0.2;

// Further constant definitions omitted ...

bl: [0..2] init 0;

el: [0..2] init 0;

[] (bl = 0) & (el = 0) -> FimB*tob: (bl’ = 1);

[] (bl = 1) & (el = 0) -> bto: (bl’ = 0) + FimB*btobb: (bl’ = 2);

[] (bl = 2) & (el = 0) -> bbtob: (bl’ = 1);

[] (el = 0) & (bl = 0) -> FimE*toe: (el’ = 1);

[] (el = 1) & (bl = 0) -> eto: (el’ = 0) + FimE*etoee: (el’ = 2);

[] (el = 2) & (bl = 0) -> eetoe: (el’ = 1);

endmodule

module rightBS

// Details omitted because of similarity to leftBS

// Essentially br replaces bl and er replaces el ...

endmodule

module switch1

// We say that state 0 corresponds to off

switch: [0..1] init 0;

[switch1] switch = 0 & ((bl = 2) & (br = 2)) ->

bbswtichon: (switch’ = 1);

[switch2] switch = 0 & ((el = 2) & (er = 2)) ->

eeswtichon: (switch’ = 1);

[switch3] switch = 1 & ((bl = 2) & (br = 2)) ->

bbswtichoff: (switch’ = 0);

[switch4] switch = 1 & ((el = 2) & (er = 2)) ->

eeswtichoff: (switch’ = 0);

endmodule

Code 6.3 A PRISM model of Fim binding

A salient feature of the graph is that, for increasing FimB, the steady state prob-

ability rapidly approaches 0.5. This can be easily explained by taking into account

that the FimB mediated switching becomes dominant as the concentration of FimB

increases. For high levels of FimB, the importance of FimE for switching becomes

negligible. That this is indeed true can be checked using transition rewards. Fig-

ure 6.7 shows the proportion of switching events that are mediated by FimB, i.e.,

that happened when all four binding sites were occupied by FimB. Clearly, as we

268 6 Other Stochastic Methods and Prism

Fig. 6.6 The steady state probability of ﬁmS being on for various concentrations of FimB

Fig. 6.7 The proportion of FimB mediated switching events in steady state

increase the concentration of FimB this proportion goes to 1. Since FimB switches

in both directions with equal rate, in the long run we would expect that the probabil-

ity of a cell to be on or off approaches 1/2 as the number of FimB increases, which

is exactly what we observe in Fig. 6.6.

6.6.2 Stochastic Versions of a Differential Equation

Differential equation models are normally much more convenient to handle and

explore than stochastic models. The formalism of differential calculus is well de-

veloped. We have also seen that it is often possible to derive exact equations for the

steady state behavior of a system, even if the full set of differential equations cannot

be solved. Unfortunately, in biology stochastic effects can often not be ignored. In

6.6 Examples 269

these situations a dual approach is helpful where a differential equation model of the

system is analyzed ﬁrst. This provides basic insights into the behavior of the system,

but tells us nothing about stochastic aspects of the system. PRISM can then be used

to simulate a stochastic version of the differential equation system.

In this second case study we give a short example of how this can be done.

Assume we have a system of differential equations for x, y, z. These variables could

be the number of particles of speciﬁc proteins in the cell, for example. We could

solve the system numerically, but that would only give us the deterministic behavior.

In order to analyze the stochastic behavior of the system we need to translate the

system of differential equations into the PRISM speciﬁcation language, which is

relatively straightforward to do. Note, however, that this will only work when the

variables of the differential equations represent numbers of entities, such as numbers

of molecules or numbers of agents. This will not always be the case. A case in point

is the differential equation model of Malaria in Sect. 4.4. If the differential equation

model represents, say, fractions of quantities, then for obvious reasons the model

cannot directly be translated into a stochastic PRISM model. Assuming this caveat

is met, we can create PRISM models directly from the differential equation (or the

differential equations) by following these steps:

1. Create a new stochastic state variable, ¯x, for each deterministic variable x in the

system of differential equations.

2. For each term in the right hand side of the differential equations, enter a new

transition in the Markov chain with the same rate as formulated in the differential

equation. For example, if the left hand side of the equation is ˙x and the right hand

side contains a term −2azx, then add the following transition:

[](x>0):2*a*z*x->(x’=x-1)

3. In the case of differential equations with several variables, there is often a ﬂow

from one variable to another, in the sense that one variable is converted into

another one. If this is the case, then this must be taken into account. For example,

if the right hand side of the equation ˙z contains a similar term +2azx, then

this term must also be taken into account in the PRISM model. The previous

transition must be amended as follows:

[] (x > 0) & (z < maxz): 2*a*z*x->

(x’=x-1)&(z’=z+1)

Note here that it was necessary to add an additional guard to prevent z from going

over its allowed limit.

Let us now illustrate this using a speciﬁc example, using a simple production-

decay differential equation.

˙x =a −bx (4.19)

270 6 Other Stochastic Methods and Prism

ctmc

const double a = 5;

const int maxx = 4000;

const double b = 0.01;

module drift

x: [0 .. maxx] init 0;

[] (x < maxx) -> a: (x’=x+1);

[](x>0)->x*b:(x’=x-1);

endmodule

Code 6.4 A production-decay model

In Chap. 4 we have shown how to derive the general solution for this equation.

Assuming initial conditions of x(0) =0 this works out to be:

x(t) =

(

1 −exp(−bt)

)

(6.54)

We do not need this general solution to generate our PRISM model, but the solu-

tion makes it more convenient to compare the behavior of the differential equation

with the model we create. (Instead, we could of course also solve (4.19) numerically,

although this seems rather silly if we have the general solution.)

Following our rules of conversion, in our PRISM version of (4.19) we need only

a single stochastic state variable, which we call x. The differential equation has a

production term and a decay term. The former translates into a transition to increase

the value of x by1witharateofa. Similarly, there is a corresponding term for

decreasing x witharateofbx. Altogether we obtain the PRISM model shown in

Code 6.4.

Figure 6.8 shows a comparison of the differential equation model with the

stochastic PRISM model. The graph illustrates well how the stochastic system ﬂuc-

tuates around the solution of the deterministic system.

6.38 Make a stochastic version of the bistable switch system (4.17) on page 178.

Compare the behavior of the stochastic and the deterministic system.

6.6.3 Tricks for PRISM Models

PRISM is a very powerful tool that is easy to use for creating and exploring stochas-

tic models in biology. Like any other software tool it does have its idiosyncrasies

6.6 Examples 271

Fig. 6.8 A stochastic simulation in PRISM compared with the differential equation

that can appear somewhat unexpected to a novice user, but they are simple to work

around if one is aware of them. We will address some of those in this section.

A common problem that can be quite puzzling at ﬁrst could be caused by a tran-

sition like this one:

[] (x < maxx) -> x * b: (x’ = x + 1)

Syntactically this line is correct. However, if a line like this is incorporated into a

model PRISM will complain about it during model building, mentioning that transi-

tion rates sum to 0. It would not build the model until the problem is ﬁxed. The root

of the problem is that the state variable x appears in the transition rate. If x =0 then

the rate would also be zero. In order to avoid this error, the guard (the conditions)

must be strengthened.

[] (x< maxx) & (x>0)-> x*b: (x’=x+1)

This version will allow the model to be built. This type of error is easy to make, but

often surprisingly difﬁcult to spot. When writing a model it is helpful to check con-

ditions carefully and to make sure that the probability cannot become zero. Some-

times, there are also situations where the transitions in the model are such that a

state variable cannot become zero; in this case PRISM may still complain if the rate

formula would formally allow negative rates.

A major limitation of PRISM when model checking is the size of the model.

There are no hard and fast rules as to what can and cannot be model checked suc-

cessfully. Generally, the more states there are, the longer the model checking takes.

272 6 Other Stochastic Methods and Prism

The speed (and feasibility) of model checking can also be inﬂuenced by the rates.

During the testing phase it is therefore important to keep the number of states as

low as possible. This allows one to obtain a feel for feasible model sizes. This

is speciﬁcally the case when using the graphical interface, which can lock during

model building or model checking, preventing the user from gracefully exiting the

program.

Model checking times can sometimes be reduced by choosing the correct op-

tions. The GUI version allows the user to select various model checking algorithms

interactively. The command line offers a number of switches to adjust how PRISM

works internally. The options can be seen by typing prism -help. A detailed

explanation of the various algorithms and options available in PRISM is beyond

the scope of this book but a summary of basic operator use may be found in Ap-

pendix A.4. Further good and practical advice is available from the FAQ section on

the PRISM web page.