Cell Cycle Control in Eukaryotes

(Lecca and Priami)

Contents

• Introduction
• The Model
• Model Checking

Introduction

In this case study we consider a probabilistic model of cell cycle control in eukaryotes, a very commonly occurring class of single-celled or multi-cellular organisms. The model is a translation of the formal specification in the biochemical stochastic π-calculus given in [LP03] which studies the relative concentration of a number types of proteins, partaking concurrently in several complex chemical reactions.

In this approach, the time until a reaction between proteins is modelled as an exponential distribution, and therefore we construct a CTMC model of the system where the states of the CTMC correspond to the number of proteins of each type and the transitions correspond to the possible reactions between the proteins. The rate of a reaction is determined by a base rate and the concentrations of the reactants (i.e. the number of each type of molecule that takes part in the reaction).

Note that, for each reaction we consider, the rate of the reaction is the base rate multiplied by the product of the number of elements of each type that take part in the reaction. Therefore, to model this is PRISM, we use the fact that in PRISM the rate of synchronous transition is defined as the product of the rates of the local transitions which form this synchronous transition.

The Model

The PRISM source code for this pathway is given below.

// cell cycle control in eukaryotes
// based on the stochastic pi caclulus model of Lecca and Priami

ctmc

// initial number of molecules
const int CYCLIN=2*N;
const int CDK=N;
const int CDH1=N;
const int CDC14=2*N;
const int CKI=N;

const int N; 

// base rates of reactions
const double R1 =0.005;
const double R2 =0.001;
const double R3 =0.003;
const double R4 =0.5;
const double R5 =0.3;
const double R6 =0.005;
const double R7 =0.009;
const double R8 =0.009;
const double R9 =0.01;
const double R10=0.017;
const double R11=0.02;

// module for including the base rates
module base_rates

	[degp]      true -> R1  : true;
	[degc]      true -> R2  : true;
	[degd]      true -> R3  : true;
	[lb]        true -> R4  : true;
	[bb]        true -> R5  : true;
	[cdh1r]     true -> R6  : true;
	[pcdh1r]    true -> R7  : true;
	[removep]   true -> R8  : true;
	[removecki] true -> R9  : true;
	[donothing] true -> R10 : true;
	[bind]      true -> R11 : true;
	
endmodule

// CYCLIN
module cyclin
	
	cyclin : [0..CYCLIN] init CYCLIN;
	cyclin_bound : [0..CYCLIN] init 0;
	degc : [0..CYCLIN] init 0;
	trim : [0..CYCLIN] init 0;
	dim : [0..CYCLIN] init 0;
	
	[lb]  cyclin>0 & cyclin_bound<CYCLIN -> cyclin : (cyclin_bound'=cyclin_bound+1) & (cyclin'=cyclin-1);
	
	[degp]  cyclin_bound>0 & degc<CYCLIN -> cyclin_bound : (degc'=degc+1) & (cyclin_bound'=cyclin_bound-1);
	[degd]  cyclin_bound>0               -> cyclin_bound : (cyclin_bound'=cyclin_bound);
	[bind]  cyclin_bound>0 & trim<CYCLIN -> cyclin_bound : (trim'=trim+1) & (cyclin_bound'=cyclin_bound-1);
	
	[degc] degc>0 -> degc : (degc'=degc-1);
	
	// rate dependent on the number of cyclin/cdk pairs which have a private channel
	[bb] trim>0 & dim<CYCLIN -> 1 : (dim'=dim+1) & (trim'=trim-1);
	
	[removecki] dim>0 & cyclin_bound<CYCLIN -> dim : (cyclin_bound'=cyclin_bound+1) & (dim'=dim-1);
	[donothing] dim>0                       -> dim : (dim'=dim);
	
endmodule

// keeps track of the number of cyclin/cdk bounded pairs (i.e. pairs what have a private channel)
module counter
	
	bound1 : [0..min(CDK,CYCLIN)]; // pairs bound (cyclin has not bound with cki)
	bound2 : [0..min(CDK,CYCLIN)]; // pairs bound (cyclin bound with cki)
	
	[lb] bound1<min(CDK,CYCLIN) -> (bound1'=bound1+1);
	[degp] cyclin_bound>0 & bound1<=cyclin_bound -> 
		bound1/cyclin_bound : (bound1'=bound1-1) // probability cyclin which is bound takes part in the reaction
		+ 1-bound1/cyclin_bound : true; // probability cyclin which is not bound takes part in the reaction
	[bind] cyclin_bound>0 & bound1<=cyclin_bound & bound2<min(CDK,CYCLIN) ->
		bound1/cyclin_bound : (bound1'=bound1-1) & (bound2'=bound2+1) // probability cyclin which is bound takes part in the reaction
		+ 1-bound1/cyclin_bound : true; // probability cyclin which is not bound takes part in the reaction
	[degc] cdk_cat>0 & bound1+bound2<=cdk_cat -> 
		bound1/cdk_cat : (bound1'=bound1-1) // probability cdk which is bound (and corresponding cyclin not bound to cki) takes part in the reaction
		+ bound2/cdk_cat : (bound2'=bound2-1) // probability cdk which is bound (and corresponding cyclin bound to cki) takes part in the reaction
		+ 1-(bound1+bound2)/cdk_cat : true; // probability cdk which is not bound takes part in the reaction
	[bb] bound2>0 -> bound2 : (bound2'=bound2-1);
	
endmodule

// CDK
module cdk
	
	cdk : [0..CDK] init CDK;
	cdk_cat : [0..CDK] init 0;
	
	[lb] cdk>0 & cdk_cat<CDK -> cdk : (cdk_cat'=cdk_cat+1) & (cdk'=cdk-1);
	[cdh1r] cdk_cat>0           -> cdk_cat : (cdk_cat'=cdk_cat);
	[degc]  cdk_cat>0 & cdk<CDK -> cdk_cat : (cdk'=cdk+1) & (cdk_cat'=cdk_cat-1);
	[bb] cdk_cat>0  -> 1 : (cdk_cat'=cdk_cat-1); // rate dependent on the number of cyclin/cdk pairs which have a private channel
	[removecki] cdk<CDK  -> 1 : (cdk'=cdk+1);
	
endmodule

// CDH1
module cdh1
	
	cdh1  : [0..CDH1] init CDH1;
	inact : [0..CDH1] init 0;
	
	[degp]    cdh1>0              -> cdh1 : (cdh1'=cdh1);
	[cdh1r]   cdh1>0 & inact<CDH1 -> cdh1 : (cdh1'=cdh1-1) & (inact'=inact+1);
	[removep] cdh1>0              -> cdh1 : (cdh1'=cdh1);
	[pcdh1r] inact>0 & cdh1<CDH1 -> inact : (inact'=inact-1) & (cdh1'=cdh1+1);
	
endmodule

// CDC14
module cdc14
	
	cdc14   : [0..CDC14] init CDC14;
	
	[pcdh1r] cdc14>0 -> cdc14 : (cdc14'=cdc14-1);
	[removep] cdc14<CDC14 -> (CDC14-cdc14) : (cdc14'=cdc14+1);
	
endmodule

// CKI
module cki
	
	cki   : [0..CKI] init CKI;
	
	[degd] cki>0 -> cki : (cki'=cki-1);
	[bind] cki>0 -> cki : (cki'=cki-1);
	
endmodule

// reward structure: abound of cdc14
rewards "cdc14"
	true : cdc14;
endrewards

// reward structure: abound of cdh1
rewards "cdh1"
	true : cdh1;
endrewards

// reward structure: abound of cyclin bound
rewards "cyclin_bound"
	true : cyclin_bound;
endrewards

The table below shows statistics for the CTMCs we have built for different numbers values of the constants N1 and N2. The tables include:

• the number of states and transitions in the CTMC representing the model;
• both the number of nodes and leaves of the MTBDD which represents the model;
• the construction time which equals the time taken for the system description to be parsed and converted to the MTBDD representing it, and the time to perform reachability, identify the non-reachable states and filter the MTBDD accordingly;
• the number of iterations required to find the reachable states (which is performed via a BDD based fixpoint algorithm).

N:		Model:			MTBDD:			Construction:
		States:	Transitions:		Nodes:	Leaves:		Time (s):	Iter.s:
2		4,666	18,342		6,880	31		0.317	27
3		57,667	305,502		22,406	48		1.627	40
4		431,101	2,742,012		81,145	78		10.025	53
5		2,326,666	16,778,785		179,015	105		32.997	66
6		9,960,861	78,768,799		352,833	139		86.077	79

Model Checking

We analyse this pathway through the following properties:

• The probability that there are k cyclin_bound elements at time T.
• The expected quantity of CDC14, CDH1 and cyclin_bound at time T

These property can be specified by the following CSL formula:

const double t; // time bound
const int k; // number of molecules

// amount of cyclin bound at time t
P=? [ true U[t,t] cyclin_bound=k ]

// expected reward at time t
R{"cdc14"}=? [I=t]
R{"cdh1"}=? [I=t]
R{"cyclin_bound"}=? [I=t]

The following graph plots, in the cases when N equals 2, 3 and 4, the probability that there are k cyclin_bound elements at time as the value of T varies.

• N = 2

  
  

• N = 3

  
  

• N = 4

  
  

In the graph below we have plotted the above expected values as T varies.

• N = 2

  

• N = 3

  

• N = 4

  

• N = 5

  

Site hosted at the Department of Computer Science, University of Oxford