https://git.exeter.ac.uk/mv286/kndy-parameter-inference
Tip revision: f5585f5f09d62a2569d21def31dc779774bc2476 authored by mv286 on 15 July 2021, 12:46:14 UTC
Deleting README.md
Deleting README.md
Tip revision: f5585f5
model_driver.m
% Function that takes care of running the model and producting the model
% prediction that will be compared with the data.
% params - a structure with the model parameters to be used
% returns retVal - a vector with the model prediction. Make sure that the
% returned vector has the same dimensions as the data
function retVal = model_driver( params)
opt = odeset('MaxStep',1);
Tspan = params.Tspan;
%% run with diestrus parameters
X0 = model_IC(params,1);
[T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,params), Tspan , X0, opt);
X0 = Y(end,:);
%calculate period and duty cycle
ft = (round(size(Y, 1)/2)):size(Y, 1);
timecourse_tmp = mean(Y(ft,3:3:end),2) ;
ac = xcorr(timecourse_tmp,timecourse_tmp);
[~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
DT = Tspan(2)-Tspan(1);
if isempty(max(diff(loc)*DT))
retVal(1, 1) = 0 ;
retVal(2,1) = 1;
else
retVal(1,1) = max(diff(loc)*DT) ;
retVal(2,1) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
end
traj(1).T = T;
traj(1).Y = Y;
%% run with estrus parameters
paramstmp = params;
paramstmp.fitparams(1:4) = paramstmp.fitparams(1:4) + paramstmp.fitparams(5:8);
X0tmp = model_IC(params,1);
[T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,paramstmp), Tspan , X0tmp, opt);
X0tmp = Y(end,:);
ft = (round(size(Y, 1)/2)):size(Y, 1);
timecourse_tmp = mean(Y(ft,3:3:end),2) ;
ac = xcorr(timecourse_tmp,timecourse_tmp) ;
[~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
if isempty(max(diff(loc)*DT))
retVal(1,2) = 0 ;
retVal(2,2) = 1;
else
retVal(1,2) = max(diff(loc)*DT) ;
retVal(2,2) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
end
traj(2).T = T;
traj(2).Y = Y;
%% run diestrus parameters + 5Hz extenal stimulation
params.p3_external = 5;
[T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,params), Tspan , X0, opt);
%calculate period and duty cycle
ft = (round(size(Y, 1)/2)):size(Y, 1);
timecourse_tmp = mean(Y(ft,3:3:end),2) ;
ac = xcorr(timecourse_tmp,timecourse_tmp);
[~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
DT = Tspan(2)-Tspan(1);
if isempty(max(diff(loc)*DT))
retVal(1,3) = 0 ;
retVal(2,3) = 1;
else
retVal(1,3) = max(diff(loc)*DT) ;
retVal(2,3) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
end
traj(3).T = T;
traj(3).Y = Y;
%% run estrus parameters + 5Hz external stimulations
paramstmp.p3_external = 5;
[T, Y] = ode23s(@(t,y) modelPopulation_v103(t,y,paramstmp), Tspan , X0tmp, opt);
%calculate period and duty cycle
ft = (round(size(Y, 1)/2)):size(Y, 1);
timecourse_tmp = mean(Y(ft,3:3:end),2) ;
ac = xcorr(timecourse_tmp,timecourse_tmp);
[~,loc] = findpeaks(ac./max(ac), 'MinPeakHeight',0.1);
DT = Tspan(2)-Tspan(1);
if isempty(max(diff(loc)*DT))
retVal(1,4) = 0 ;
retVal(2,4) = 1;
else
retVal(1,4) = max(diff(loc)*DT) ;
retVal(2,4) = mean(timecourse_tmp > ((max(timecourse_tmp) -min(timecourse_tmp))/2));
end
traj(4).T = T;
traj(4).Y = Y;
end
%ydata = [1/19.257];% 0.149]