%In this script, we will set up a numerical integration of a genetic
%oscillator. It's important to note that the oscillation only occurrs under
%a certain set of conditions which are not inherently obvious. You should
%try playing with these parameters (decay rate, generation rate, and
%initial conditions) to see what permits and what breaks the oscillatory
%behavior. As always, we'll start by setting a few variables.
r0A = 100; %Basal rate of activator A expression (1/min).
r0R = 1; %Basal rate of repressor R expression (1/min).
rR = 100; %Expression of R when activated by A (1/min).
rA = 5000; %Expression of A when activated by A (1/min).
gammaA = 30; %Decay rate of A (1/min).
%To set up the numerical integration, we need to figure out how long to run
%the simulation and what timestep we will take.
TotalTime = 20; %Total length of the simulation (in min).
nSteps = 1000; %Number of steps to integrate.
dt = TotalTime / nSteps; %Time steps.
%We should set up empty vectors to which we'll add the amount of A and R at
%each time step.
cA = zeros(nSteps,1);
cR = zeros(nSteps,1);
time = zeros(nSteps,1); %Also to keep track of how much time has elapsed.
%The oscillatory behavior of this genetic network is very dependent on the
%starting conditions. This will oscillate with cA = 1 and cR = 10. You
%should try to find conditions where you lose the oscilliatory behavior.
cA(1) = 1; %Initial concentration of A activator (1/cell).
cR(1) = 10; %Intial concentration of R repressor (1/cell).
%Now we can set up the numerical integration. This is basically the exact
%same thing we did when we integrated the ligand-receptor binding.
for i=1:nSteps - 1
%Let's compute the change in cA and cR at each i'th time point.
dcA = -gammaA*cA(i) + (r0A + rA*cA(i)^2) / (1 + cA(i)^2 + cR(i)^2);
dcR = -cR(i) + (r0R + rR*cA(i)^2) / (1 + cA(i)^2);
%Now we'll use this to figure out the concentration of repressor and
%activator at the i + 1 time point.
cA(i+1) = cA(i) + dcA*dt;
cR(i+1) = cR(i) + dcR*dt;
%We can keep track of time by adding the time step at each point.
time(i+1) = time(i) + dt;
end
%Now, let's see if it oscillates.
plot(time, cA, 'r-')
hold on
plot(time, cR, 'b-')
legend('Activator A', 'Repressor R', 'Location', 'northeast');
xlabel('Time (min)');
ylabel('Protein Concentration (# per cell)');
hold off