function [t,y] = MOtesterAP_post(m, thor, endTime)
%
% Dynamic simulation that uses DDE to follow her/hes mRNA and protein in a single cell over time
% in first the posterior, than anterior PSM. Requires three inputs: m (which knockdown condition)
% thor (GRN) and endTime (total time for simulation, usually 300 minutes)
%
% Michael Sneddon 2007
% Jamie Schwendinger-Schreck 2013
%%% Knockdown conditions for all m
% m=1 her1 kd
% m=2 her7 kd
% m=3 her11 kd
% m=4 her12 kd
% m=5 hes6 kd
% m=6 her15 kd
% m=7 her1+7 kd
% m=8 her1+11 kd
% m=9 her1+12 kd
% m=10 her1+6 kd
% m=11 her1+15 kd
% m=12 her7+11 kd
% m=13 her7+12 kd
% m=14 her7+6 kd
% m=15 her7+15 kd
% m=16 her11+12 kd
% m=17 her11+6 kd
% m=18 her11+15 kd
% m=19 her12+6 kd
% m=20 her12+15 kd
% m=21 hes6+her15 kd
% m=22 wild type
% Parameters for the simulation
startTime = 0;
delayStep = .1;
initValue = 0;
% Define a time array and a value array (to store past values of a variable
% For multiple variables, valueArray should be a matrix.
timeArray = [startTime:delayStep:endTime]';
valueArray = zeros(13, length(timeArray))+initValue;
curIndex = 1;
% Parameters for the model
% Define constants
deltaM = 0.31; %deg rate of her mRNA
deltaM6 = 0.31; %deg rate of hes6 mRNA
deltaP = 0.15; %deg rate of Her protein
deltaP6 = 0.15; %deg rate of Hes6 protein
deltaD = 0.15; %deg rate of Her1 dimers
deltaD6 = 0.15; %deg rate of Hes6/Hes6 dimer
betaM = 150; %max rate of synthesis for her mRNA
betaM6 = 150; %max rate of synthesis for hes6 mRNA
%betaP=matrix(22,6);
%where m calls row of betaP values, corresponding to some her/hes knockdown
betaP =[ 0.017,1,1,1,1,1;
1,0.015,1,1,1,1;
1,1,0.095,1,1,1;
1,1,1,0.019,1,1;
1,1,1,1,0.03,1;
1,1,1,1,1,0.022;
0.017,0.015,1,1,1,1;
0.017,1,0.095,1,1,1;
0.017,1,1,0.019,1,1;
0.017,1,1,1,0.030,1;
0.017,1,1,1,1,0.022;
1,0.015,0.095,1,1,1;
1,0.015,1,0.019,1,1;
1,0.015,1,1,0.03,1;
1,0.015,1,1,1,0.022;
1,1,0.095,0.019,1,1;
1,1,0.095,1,0.03,1;
1,1,0.095,1,1,0.022;
1,1,1,0.019,0.03,1;
1,1,1,0.019,1,0.022;
1,1,1,1,0.03,0.022;
1,1,1,1,1,1];
Tns = 18.5; %translation rate of Her proteins
Tns6 = 18.5; %translation rate of Hes6
betaP1 = Tns*betaP(m,1);
betaP7 = Tns*betaP(m,2);
betaP11 = Tns*betaP(m,3);
betaP12 = Tns*betaP(m,4);
betaP6 = Tns6*betaP(m,5);
betaP15 = Tns*betaP(m,6);
%where m=22 indicates wt values, m=1 indicates her1 mo inj, and so on...
r = 13900; %critical level of Her dimer for repression of her txn
%equals sum of all r values for Cinquin, divided by 7
%(#dimer represented) Same for all genes.
s = 162.7; %level of repression associated with Her1-1 dimer on her1 txn
%repression is further scaled by "thor" values, where a
%stron dimer is s*2, a weak dimer is s and no repression is
%s*0.
% s calculated by summing all s values in
%Cinquin, 2007; it should then be divided by total possible number
%of repressing dimers (i.e. in post vs ant).
kOn = 0.02; %association rate of Her dimers
kOff = 1.0; %dissociation rate of Her dimers
kOn6 = 0.02; %association rate of Hes6/Hes6
kOff6 = 1.0; %dissociation rate of Hes6/Hes6
delaym = 5.0; %delay associated with txn, splicing, export
delay = 1.5; %delay associated with translation, import to nucleus
% Starting concentrations
initConc = zeros(22,1);
% SEDIT - made an array with the list of her's that is in the same order as
% the function xdot.
herNames = ['her1 mRNA','her7 mRNA','her11 mRNA','her12 mRNA','Hes6 mRNA', ...
'her15 mRNA','Her1 protein','Her7 protein','Her11 protein', ...
'Her12 protein','Hes6 protein','Her15 protein','Her1-1 Dimer', ...
'Her7-6 Dimer','Her11-11 Dimer','Her12-12 Dimer','Her12-6 Dimer', ...
'Hes6-6 Dimer (non DNA-binding)','Her6-15 Dimer','Her15-15 Dimer', ...
'Her1-6 Dimer (non DNA-binding)','Her7-15 Dimer (non DNA-binding)'];
% Simulate the system as you normally would
[t,y] = ode15s(@xdot,timeArray,initConc);
% Function for assessing the past values [Michael Sneddon, 2007]
function [pastVal] = getPastValue(t)
% first find the index of the time point given, which
% will give you the index=1 if you picked a time earlier
% than any point in the timeArray...
ind=find(timeArray>=t,1,'first');
% Check that we found a value (if we did not, then we are trying
% to lookup a time that exceeds the timeArray
if isempty(ind)
fprintf('ERROR: in getPastValue(t), t exceeds time array.\n');
pastVal=0;
return;
end
% First approximation: we can just return the value
% directly at that time
pastVal = valueArray(:, ind);
% Better approximation: interpolate between the points
% given to get a better estimate of what the value should be
%
% pastVal = interp( ... )
end
function [df]=xdot(t,f)
%fprintf(['t=',num2str(t),'\n']);
h1m = f(1);
h7m = f(2);
h11m = f(3);
h12m = f(4);
h6m = f(5);
h15m = f(6);
h1 = f(7);
h7 = f(8);
h11 = f(9);
h12 = f(10);
h6 = f(11);
h15 = f(12);
h1_1 = f(13);
h7_6 = f(14);
h11_11 = f(15);
h12_12 = f(16);
h12_6 = f(17);
h6_6 = f(18);
h6_15 = f(19);
h15_15 = f(20);
h1_6 = f(21);
h7_15 = f(22);
% be sure to set our array in here, say if the value we want
% to keep track of is at position one
ind=find(timeArray>=t,1,'first');
if(curIndex