This example illustrate the use of PRISM for a small biological case study. We model circadian rhythms which are used by a wide range of organisms, from bacteria to mammals, to keep a sense of daily time and regulate behaviour accordingly. In particular, we use the abstract circadian clock model of Barkai and Leiber. See the following links for background reading on this example: [BL00, VKBL02].
The system comprises seven different species:
The following summarises the reactions that occur in the system:
A PRISM model of this system is shown below. This is a continuous-time Markov chain (CTMC), which represents a discrete stochastic model of the system. For further reading on this topic, see for example section 2 (and in particular section 2.2) of [KNP+06].
// stochastic model of circadian clock based on ODE model of // J. Vilar et al. Mechanisms of Noise-Resistance in Genetic Oscillators // Proc. National Academy of Sciences of the United States of America, 99(9):5988-5992, 2002 // dxp/gxn 24/04/08 ctmc // generic bound (for species with potentially unbounded quantities) const int N=1000000; //-------------------------------------------------------- // activator modules module activator_gene da : [0..1] init 1; // amount of activator gene da_a : [0..1]; // amount of activator gene with activator protein bound [transc_da] da=1 -> da : true; // transcription (gene to mRNA) [transc_da_a] da_a=1 -> da_a : true; // transcription (gene to mRNA) (activator protein bound to gene) [bind_a] da=1 & da_a=0 -> da : (da'=0) & (da_a'=1); // binds with activator protein [rel_a] da=0 & da_a=1 -> da_a : (da'=1) & (da_a'=0); // unbinds from activator protein endmodule module activator_mRNA ma : [0..N]; // amount of activator mRNA [transc_da] ma<N -> (ma'=ma+1); // transcription (gene to mRNA) [transc_da_a] ma<N -> (ma'=ma+1); // transcription (gene to mRNA) (activator protein bound to gene) [transl_a] ma>0 -> ma : true; // translation (mRNA to protein) [deg_ma] ma>0 -> ma : (ma'=ma-1); // mRNA degrades endmodule module activator_protein a : [0..N]; // amount of activator protein [transl_a] a<N -> (a'=a+1); // translation (mRNA to protein) [bind_a] a>0 -> a : (a'=a-1); // protein binds with activator gene [bind_r] a>0 -> a : (a'=a-1); // protein binds with repressor gene [rel_a] a<N -> (a'=a+1); // unbinds from activator gene [rel_r] a<N -> (a'=a+1); // unbinds from repressor gene [deg_a] a>0 -> a : (a'=a-1); // protein degrades [deactive] a>0 -> a : (a'=a-1); // deactivation (binds with repressor protein) endmodule //-------------------------------------------------------- // repressor modules // rename activator gene module module repressor_gene = activator_gene[da=dr, da_a=dr_a, bind_a=bind_r, rel_a=rel_r, transc_da=transc_dr, transc_da_a=transc_dr_a] endmodule // rename activator mRNA module module repressor_mRNA = activator_mRNA[ma=mr, transl_a=transl_r, transc_da=transc_dr, transc_da_a=transc_dr_a, deg_ma=deg_mr] endmodule module repressor_protein r : [0..N]; // amount of repressor protein [transl_r] r<N -> (r'=r+1); // translation (mRNA to protein) [deg_r] r>0 -> r : (r'=r-1); // protein degrades [deactive] r>0 -> r : (r'=r-1); // deactivation of activator protein (protein binds with activator protein) [deg_c] r<N -> (r'=r+1); // activator protein degrades (when repressor protein is bound) endmodule //-------------------------------------------------------- // inactive protein module module inactive_protein c : [0..N]; // amount of inactivated protein [deactive] c<N -> (c'=c+1); // protein forms (activator and repressor proteins bind) [deg_c] c>0 -> c : (c'=c-1); // activator protein degrades endmodule //-------------------------------------------------------- // stochastic rates // rates module module rates [transc_da] true -> 50 : true; // transcription rate of activator gene [transc_da_a] true -> 500 : true; // transcription rate of activator gene (activator protein bound) [transc_dr] true -> 0.01 : true; // transcription rate of repressor gene [transc_dr_a] true -> 50 : true; // transcription rate of repressor gene (activator protein bound) [transl_a] true -> 50 : true; // translation rate of activator mRNA [transl_r] true -> 5 : true; // translation rate of repressor mRNA [bind_a] true -> 1 : true; // binding rate of activator gene and activator protein [bind_r] true -> 1 : true; // binding rate of repressor gene and activator protein [deactive] true -> 2 : true; // binding rate of activator and repressor proteins [rel_a] true -> 50 : true; // release rate of activator gene and activator protein complex [rel_r] true -> 100 : true; // release rate of repressor gene and activator protein complex [deg_a] true -> 1 : true; // degradation rate of activator protein [deg_c] true -> 1 : true; // degradation rate of activator protein [deg_r] true -> 0.2 : true; // degradation rate of repressor protein [deg_ma] true -> 10 : true; // degradation rate of activator mRNA [deg_mr] true -> 0.5 : true; // degradation rate of repressor mRNA endmodule // amount of activator protein rewards "activated_protein" true : a; endrewards rewards "repressor_transcriptions" [transc_dr] true : 1; [transc_dr_a] true : 1; endrewards
The PRISM module comprises 7 modules. Each one represents the current quantity of one of the 7 species of the system, as listed above. For example, module activator_protein has a single variable a, representing the number of activator proteins. The range of values for this variable is [0..N]. If you look at the top of the file, you will see that N, which is used as the upper bound for most species, has been set to 1000000. This is because, although there are potentially an unbounded number of these species, PRISM requires models to be finite state. Due to the size of this model we will only use the simulator engine during our analysis.
Look at the module activator_gene. It has a 0-1 variable da, representing the absence (0) or presence (1) of an activator gene (for now, we assume that only one such gene can be present). This module has an additional variable da_a, which we use to denote the fact that the gene is bound to an activator protein (A).
The modules representing genes and mRNA for the repressor are identical in structure to those for the activator and so are created with module renaming. The initial state of the model is given by the initial state of each variable. For da (and thus also dr), this is 1. For all other variables, the initial value is taken by default to be the lowest value in its range. Thus, in the initial configuration of our model, there is just a single activator gene and a single repressor gene.
Download the model file circadian.sm from above and load it into PRISM.
Create a new path in the simulator and look at the values of the variables
in the initial state. These should correspond to the description above.
Look at the "Manual exploration" box at the top
to see what transitions are available from the initial state.
We will now see how the first of these, labelled transc_da,
relates to the model description.
When reactions occur in the system, the variables representing the amount of each species will be updated accordingly. We use synchronisation to allow the modules representing different species to update simultaneously. Look at the two commands in the model file above that are labelled with action transc_da: one in module activator_gene and one in activator_mRNA (action labels appear in the pair of square brackets at the start of a command).
For a transc_da-labelled transition to occur, the guards for both of these commands must be satisfied: da=1 (the activator gene is present) and ma<N (the amount of activator mRNA has not reached its maximum). When the transition occurs, the right-hand sides of both commands are applied. For the first command, the right-hand side is just true, indicating that no change occurs. For the second, the amount of activator mRNA is incremented: (ma'=ma+1).
The rate of this synchronous transition is determined by multiplying the rates of each command that contributes to it. Rates are assigned to commands in the same way that probablities are. For example, the 3rd and 4th commands of module activator_mRNA have rate ma. In commands where the rate is omitted, it is assumed to be 1. For the transc_da-labelled transition we have just described, both commands have rate 1. So where does the rate of 50.0 (that you can see in the simulator) come from? If you look at the end of the model, you will see a module called rates containing a number of "dummy" commands whose only purpose is to contribute to the definition of the rate of synchronous transitions. Hence, transc_da transitions have rate equal to 50.
Use the PRISM simulator to step through some random paths through the model.
For each step, try to understand its corresponding definition in the PRISM model,
and also what it relates to in the high-level description of the circadian clock,
given at the top of this page.
Generate a path in the simulator where the number of mRNA of A reaches 5.
What happens to the rate of translation of A
(action transl_a) as the number of mRNA increases?
We will now use PRISM to analyse the behaviour of the circadian clock. The crucial feature of circadian clocks is that the clock runs accurately, triggering the appropriate genes at the correct times of the day. This feature is ensured through the clock's ability to maintain a constant period over a wide range of external and internal fluctuations.
More specifically, we will first study the quantity of the activator proteins as time evolves. To do this we add the following reward structure to the PRISM description.
// amount of activator protein rewards "activated_protein" true : a; endrewards
Read the section on
costs and rewards
in the manual.
Then, add this reward structure to your PRISM description of the circadian clock.
Go to the "Properties" tab of the GUI,
and add a new constant called T
in the "Constants" panel.
T should be of type double and have no defined value.
Then add the following property (the expected value of the reward structure
"activated_protein" at time instant T):
R{"activated_protein"}=? [ I=T ]
Now, create an experiment for the property,
for T from 0 up to 200,
i.e. for the first 200 hours of the system.
Important:
Do this using the simulation engine (i.e. tick the "Use simulation" box in the "Define constants" window),
set the number of samples to be generated by the simulator to 1,
and increase the "Maximum path length" setting to 1000000.
Look at the resulting graph that is generated
(you may may have to wait for a few seconds for simulation to complete).
You should observe oscillations in the amount of activator proteins,
occurring at regular time intervals.
Next, we will study the influence of the rates for transcription on the oscillation of the activator protein.
Add a constant called
k of type double to the model and change the
rates of transcription such that the rates are multiplied by a factor of k. Check your modifications are correct by generating
a number of random paths in the simulation engine.
Run an experiment, similar to the one above, to see what happens to the
oscillations of the activator protein as k
varies between 0.5 and 1.5.
Make sure that you choose T as the variable
to be plotted on the x-axis in the resulting graph.
When prompted, say that you want to use the same initial state for the experiments). The results should
demonstrate that changing the rate of transcription influences the peaks of the
oscillations but not the period. Note that the peaks will not be at the same
points due to the stochastic nature of the model and the fact the plots are
for different simulation runs.
We will now look at the number of transcriptions that occur over a given period of time.
Add a reward structure
to the model called "repressor_transcriptions",
which assigns 1 to transitions which
corresponds to a transcription of the repressor gene.
Now create a new property to check the expected cumulated
reward up until time T.
Don't forget to specify which reward structure you want in the property.
(You might want to look at
this section
of the manual.)
How does this measure vary for different values of k?
Finally, we will look at a more general model where initially there is more than one activator gene and more than one repressor gene.
Add an integer constant M to the model and rewrite the
activator_gene module to model the case
where there are initially M activator genes. Note that since the
repressor gene is constructed by renaming the activator gene you are in fact changing the
model so that the initial number of repressor genes is also M.
To ensure that you have modified the model correctly, use the simulation engine to generate some
random paths and check that the behaviour is as you expect.
What effect does increasing the initial numbers of genes to 2, 3 and 4 have
on the oscillations of the activator protein?
You can fix k at 1 for this experiment.
The results should again show that the clock runs accurately under fluctuations: changing
the initial numbers of genes does influence the values of the peaks, but does not
change the period of the oscillations.
[ Back to index ]