%Simulate mRNA birth and death using the Gillespie algorithm
%Model parameters
r=3; %1/min
gamma=1/3; %1/min/mRNA
%Simulation parameters
TotalReactions=100; %Total number of events we want to simulate
Time=0; %This initializes a vector where we will keep track of time
%Initial condition
m(1)=0;
for i=2:TotalReactions %Repeat the loop TotalReactions times
%First, get the probabilities per unit time corresponding to
%this time point
k1=r; %Production event
k2=gamma*m(i-1); %Degradation event
k0=k1+k2;
%Get a random number to decide on the DeltaT to the next
%reaction
Coin1=rand;
%Weight this uniformly distributed random number using an
%exponential distribution
DeltaT=1/k0*log(1/Coin1);
%Add this DeltaT to the time of the previous event to get the
%time of the current event
Time(i)=DeltaT+Time(i-1);
%Second, flip a new coin to decide which event will happen
Coin2=rand;
if Coin2