%% This function uses a Finite State Projection to analyse a model of the %% genetic toggle switch. This code is to be used in %% conjunction with homework assignment. This is by no means an optimized code, %% and will require some modifications to fully complete the homework. File %% created by brian munsky and made availalble to Q bio summer school %% students, in July, 2008. clear all %% Toggle Paramters alpha1=16; alpha2=50; beta = 2.5; gamma = 1; %% Choose the Configuration Space to include. N=[40,100]; %Max (kept) numbers for species 1 and species 2 %% Initiaize the generator matrix, A. lenA = (N(1)+1)*(N(2)+1); %this is the number of points in the %configuration rectangle, plus two absorbing sinks A = spalloc(lenA+2,lenA+2,0); %A is created with sparsity structure and set to zero. %% Now to populate the matrix A. Ikeep=[]; %Kepp track of states we will later want to keep. for i=0:N(1) %number of species 1 for j=0:N(2) % number of species 2 ix = i*(N(2)+1)+j+1; % index of current state in rectangle. if i*j<=260 % we are going to exclude all states such that x1x2>260 % For the later parts of the homework, this part will % need to be removed, or you can simply set it to a % much larger cut off. Ikeep=[Ikeep,ix]; %Keep the rest. else break %(optional) if the current state is outside our region of % interest, j+1 will also be. end ix_2p = ix+1; %index of x + s1 in rectangle ix_2d = ix-1; %index of x + s2 in rectangle ix_1p = ix+(N(2)+1); %index of x + s3 in rectangle ix_1d = ix-(N(2)+1); %index of x + s4 in rectangle prop = [alpha1/(1+j^gamma),... i,... alpha2/(1+i^beta),... j]; %Propensities of each reaction. % first Reaction. A(ix,ix) = A(ix,ix)-prop(1); %negative diagonal element if i0 A(ix_1d,ix) = prop(2); end %third reaction A(ix,ix) = A(ix,ix)-prop(3); if j0 A(ix_2d,ix) = prop(4); end end end Ikeep=[Ikeep,length(A)-1,length(A)]; % we also want to keep the two sink terms. A = A(Ikeep,Ikeep); %find the FSP generator with the states we chose. P0 = zeros(length(A),1); %initialize the probability distribution: s1=s2=0 with Prob 1; P0(1)=1; %% *********************************************************** % Now that we have the generator, lets plot some results! %*********************************************************** % First lets look at a probaility distribution at a specific time t. tau = 10; %time at which to compute the distribution. ExpAt = expm(A*tau); %Here I use a matrix exponential. There are much faster approaches. %% P_tau(Ikeep) = ExpAt*P0; %Compute P at time tau, and transfer back to the % rectangular coordinates (for plotting purposes) P3d_tau = spalloc(N(2)+1,N(1)+1,0); for i=0:N(1) %number of species 1 for j=0:N(2) % number of species 2 ix = i*(N(2)+1)+j+1; % index of current state in rectangle. P3d_tau(j+1,i+1) = P_tau(ix); end end XX = [0:N(1)]; YY = [0:N(2)]; figure(1) mesh(XX,YY,log(P3d_tau)) return %% Next lets look at the flow of probability exit into the absorbing states as %% functions of time. Hint: Something like this is needed for the latter parts of the %% homework. You should be able to choose the right projections (N) yourself!! t=0; % initial time tt=[]; %vector of times xx=[]; %vector of the two absorbing sink probailities. P= P0; %initial probability distribution from above. for j=-1:2 %I will use different time scales of sizes 10^j tau = 10^j; %Time step ExpAt = expm(A*tau); %FSP map for given time step. for i=1:300 P = ExpAt*P; %Update dist. t=t+tau; %update time. tt=[tt,t]; %record time xx=[xx,P(length(P)-1:length(P),1)]; %record prob in both sinks. end end figure(2) loglog(tt,xx(1,:),'g','linewidth',2) hold on loglog(tt,xx(2,:),'r','linewidth',2)