Skip to content

GitLab

Commit 067d60d4 authored by Antonio Politi's avatar Antonio Politi
Browse files

Modeling variants equations 

* Alternative production
* Multistep and distributed delay model
parent 26cede72
Hide whitespace changes
Inline Side-by-side
matlabcode/delayedEquationModel/absDelayedMaturation.m
View file @ 067d60d4
%%
% abstract Class for delayed pore maturation
%
classdef absDelayedMaturation < handle
properties (Access = protected)
timeI = [0:150]; %intergration time points
dataA; %data structure containing array of all data
tpm; %time of postmitotic start for non-core NC and Core region [time_start_pm_NC time_start_pm_C]
ti; %time interphase assembly start [time_start_ip_NC time_start_ip_C]
dataA; %data structure containing array of all data
end
properties
indir;
outdir;
properties
name; % name of the model. This is the name of the class
indir; % input directory of the data
outdir; % output directory of the data
dataF; %data structure containing data to be fitted
useSpline; % use spline fit to surface area or combination of exponential and linear function
useStd = 0; % scale each data point with its std
useIntermediate = 1; % if 0 fit only mature pore
LB; % Lower boundary for fitting data
HB; % Higher boundary for fitting data
parnames; % name of parameters
N; % number of variables
end
methods
function MO = absDelayedMaturation(force, showplot)
methods (Abstract)
[soli, solo] = solve(MO, par, timeint) %implements the solver of the model
[obs] = getObservable(MO,sol) %implements the observable of the model (intermediates and mature pores)
[prodRate, maturePM] = getProdRate(MO, par, timeint)
end
methods
%%
% absDelayedMaturation(surf, useSpline, force)
% constructor
% surf: gives which surface data should be used. 24h, or 2h
% useSpline: specify if spline fit of surface is used
% force: force reading in of input data
function MO = absDelayedMaturation(surf, useSpline, force)
if nargin < 1
force = 0;
surf = '24h';
end
if nargin < 2
showplot = 0;
MO.useSpline = 0;
else
MO.useSpline = useSpline;
end
if nargin < 3
force = 0;
end
MO.dataA = MO.dataIn(surf, force);
end
%%
% fitModel(ifit, par)
% Run fit of model
% ifit: index of parameter(s) to fit
% par: parameter values
function [parf, norm] = fitModel(MO, ifit, par)
lb = MO.LB(ifit);
hb = MO.HB(ifit);
options = optimset('Display', 'iter');
[parf, norm] = lsqnonlin(@MO.distMO, par(ifit), lb, hb, ...
options, ifit, par);
end
%%
% distMO(MO, parf, ifit, par)
% compute vector of distance model to data
% parf: value of parameters to be fitted
% ifit: index of parameters to be fitted
% par: vector containing all parameters
function [dist, distI, distM] = distMO(MO, parf, ifit, par)
parloc = par;
parloc(ifit) = parf;
[soli, solo] = solve(MO, parloc, [par(1) 120]);
A = MO.getSurf(MO.dataF(1).time);
yi = deval(soli, MO.dataF(1).time);
yo = deval(solo, MO.dataF(1).time);
obsi = MO.getObservable(yi);
obso = MO.getObservable(yo);
distI = obsi(:,1)./A - MO.dataF(1).core_i(:,1);
distI = [distI; obso(:,1)./A - MO.dataF(1).core_o(:,1)];
distM = obsi(:,2)./A - MO.dataF(2).core_i(:,1);
distM = [distM; obso(:,2)./A - MO.dataF(2).core_o(:,1)];
if MO.useIntermediate
dist = [distI;distM];
stdVec = [MO.dataF(1).core_i(:,2);MO.dataF(1).core_o(:,2);MO.dataF(2).core_i(:,2);MO.dataF(2).core_o(:,2)];
else
dist = [distM];
stdVec = [MO.dataF(2).core_i(:,2);MO.dataF(2).core_o(:,2)];
end
if MO.useStd
distI = distI./stdVec(1:length(distI));
distM = distM./stdVec(length(distI)+1:length(distM));
dist = dist./stdVec;
else
distI= distI/mean(stdVec(1:length(distI)));
if length(stdVec) < length(distI)+length(distM)
distM = distM/mean(stdVec);
else
distM = distM/mean(stdVec(length(distI)+1:end));
end
dist = dist/mean(stdVec);
end
end
%%
% runFit(par, ifit, maxNrfit, fact, outfile)
% performs multiple fits and write out the results
% par: Values of all paramaters
% ifit: indexes of parameters to fit
% maxNrfit: maximal number of rounds
% fact: set the initial values of the parameters to fit
% (randomly distributed par_ini +- fact*rand*par_ini)
% outfile: name of output file
function [par, norm] = runFit(MO, par, ifit, maxNrfit, fact, outfile)
try
parFileName = fullfile(MO.outdir, outfile); % full output file name
try
parfitT = load(parFileName);
inifit = size(parfitT,1); %where to start the fitting procedure
fid = fopen(parFileName, 'a');
catch
fid = fopen(parFileName, 'w');
fprintf(fid,'%% Model %s\n',MO.name);
fprintf(fid,'%% N %d\n', MO.N);
parnamesA = strsplit( MO.parnames, ' ');
fprintf(fid, '%% useStd %d useIntermediate %d useSpline %d\n', MO.useStd, MO.useIntermediate, MO.useSpline );
fprintf(fid, '%% fittedIds %s\n', num2str(ifit));
fprintf(fid, '%% fittedPar %s\n', strjoin(parnamesA(ifit), ' '));
fprintf(fid, '%% factParvariability %.2f\n', fact);
fprintf(fid, '%% %s\tnorm\tnrpara\n', strjoin(parnamesA,'\t'));
inifit = 0;
end
for nrfit = inifit+1:maxNrfit
parfit = par;
parfit(ifit) = parfit(ifit)+ fact*parfit(ifit).*(1-2*rand(1,length(ifit)));
parfit = (parfit > MO.HB).*(MO.HB) + (parfit <= MO.LB).*(MO.LB) + (parfit <= MO.HB).*(parfit > MO.LB).*parfit;
[parf, norm] = MO.fitModel(ifit, parfit);
parfit(ifit) = parf;
for i=1:length(parfit)-1
fprintf(fid, '%.3e\t', parfit(i));
end
fprintf(fid, '%.3e\t%.3e\t%d\n', parfit(end), norm, length(ifit));
end
fclose(fid);
catch
fclose(fid);
end
fclose all;
end
function par = getBestfit(MO, filename)
try
parV = load(filename);
catch
error(sprintf('could not load %s\n', filename));
end
if parV(1, end) == parV(2,end)
[val, ip] = min(parV(:,end-1));
par = parV(ip, 1:end-2);
else
[val, ip] = min(parV(:,end));
par = parV(ip, 1:end-1);
end
end
function Den = writeDensities(MO, par, outfile)
timeplot = [par(1):0.1:120]';
A = MO.getSurf(timeplot);
[soli, solo] = solve(MO, par,[par(1) 120]);
yi = deval(soli, timeplot);
yo = deval(solo, timeplot);
obsi = MO.getObservable(yi);
obso = MO.getObservable(yo);
Den = [timeplot obsi(:,1)./A obsi(:,2)./A obso(:,1)./A obso(:,2)./A];
fid = fopen(outfile, 'w');
if fid > 0
fprintf(fid,'%s\t%s\t%s\t%s\t%s\n','time', 'inter_i', 'mature_i', 'inter_o', 'mature_o')
for i=1:length(timeplot)
fprintf(fid, '%d\t%.2e\t%.2e\t%.2e\t%.2e\n',Den(i,1), Den(i,2),Den(i,3),Den(i,4),Den(i,5));
end
fclose(fid);
else
error('cant open outputfile');
end
MO.dataA = MO.dataIn(force, showplot);
end
function dataA = dataIn(MO, force, showplot)
function writeProdRate(MO, par, outfile, timeint)
if nargin < 4
timeint = [10:0.1:120]';
end
[prodrate, maturePM] = MO.getProdRate(par, timeint);
fid = fopen(outfile, 'w');
fprintf(fid, 'time\tProdDensity_inner\tProdDensity_outer\n');
timeint = [0;9.95;10; timeint];
prodrate = [[0;0;par(end-3)/0.1/MO.getSurf(10); prodrate(:,1)] [0;0;par(end-1)/0.1/MO.getSurf(10); prodrate(:,2)] ];
for i=1:length(timeint)
fprintf(fid, '%.2f\t%.3e\t%.3e\n',timeint(i), prodrate(i,1), prodrate(i,2) );
end
fclose(fid);
end
function writeFractionMatureIP(MO, par, outfile, timeint)
if nargin < 4
timeint = [10:0.1:1200]';
end
[soli, solo] = MO.solve(par, [par(1) timeint(end)]);
[prodRate, maturePM] = MO.getProdRate(par, timeint);
A = MO.getSurf(timeint');
yi = deval(soli, timeint);
yo = deval(solo, timeint);
maturei = (yi(3,:)./A)';
matureo = (yo(3,:)./A)';
fid = fopen(outfile, 'w');
fprintf(fid, 'time\tFractionMatureIpInner\tFractionMatureIpOuter\n');
for i=1:length(timeint)
fprintf(fid, '%.2f\t%.3e\t%.3e\n',timeint(i), (maturei(i)- maturePM(i,1))/maturei(i), ...
(matureo(i) - maturePM(i,2))/matureo(i) );
end
fclose(fid)
end
%%
function plotBestfit(MO, filename, ifig, save)
if nargin < 3
ifig = 1;
end
MO.plotResult(MO.getBestfit(filename), ifig);
if save
[pathstr, filen, exte] = fileparts(filename);
saveas(ifig, fullfile(pathstr, [filen '_inner.pdf']));
saveas(ifig+1, fullfile(pathstr, [filen '_outer.pdf']));
end
end
%%
% plot the model and data given the paramter par
function [par, norm] = plotResult(MO, par, ifig)
if nargin < 3
ifig = 1;
end
set(0,'DefaultTextInterpreter','none');
set(0,'defaultAxesFontName', 'Arial');
[soli, solo] = solve(MO, par,[par(1) 1200]);
timeplot = [par(1):1:120 ]';
A = MO.getSurf(timeplot);
yi = deval(soli, timeplot);
yo = deval(solo, timeplot);
obsi = MO.getObservable(yi);
obso = MO.getObservable(yo);
barwidth = 1;
figsize = [550 350]/2;
Mpc = [30 71 126]/255;
mpc = [237 30 36]/255;
Mpc_mod = Mpc;
mpc_mod = mpc;
%%
setupFigure(ifig, [50 600 figsize],'on', 'pixels', 8);
plot(timeplot, obsi(:,1)./A, '--', 'LineWidth', 2, 'Color', mpc_mod )
plot(timeplot, obsi(:,2)./A, '--', 'LineWidth', 2, 'Color', Mpc_mod );
hb = bar(MO.dataF(1).time, [MO.dataF(2).core_i(:,1) MO.dataF(1).core_i(:,1) ], 'BarWidth', barwidth)
xlim([0 120]);
ylim([0 13]);
set(hb(2), 'FaceColor', mpc, 'EdgeColor', mpc)
set(hb(1), 'FaceColor', Mpc, 'EdgeColor', Mpc)
xlabel('Time after AO (min)')
ylabel('Density (pores/um2)')
set(gca, 'YTick', [0 4 8 12], 'XTick', [0:20:120], 'TickDir','out', 'LineWidth', 1.5,'FontName', 'Arial')
setupFigure(ifig+1, [50 600 figsize],'on', 'pixels', 8);
outerD = [obso(:,1)./A obso(:,2)./A];
plot(timeplot, obso(:,1)./A, '--', 'LineWidth', 2, 'Color', mpc_mod )
plot(timeplot, obso(:,2)./A, '--', 'LineWidth', 2, 'Color', Mpc_mod );
hb = bar(MO.dataF(1).time, [MO.dataF(2).core_o(:,1) MO.dataF(1).core_o(:,1) ], 'BarWidth', barwidth)
xlim([0 120]);
ylim([0 13]);
set(hb(2), 'FaceColor', mpc, 'EdgeColor', mpc)
set(hb(1), 'FaceColor', Mpc, 'EdgeColor', Mpc)
xlabel('Time after AO (min)')
ylabel('Density (pores/um2)')
set(gca, 'YTick', [0 4 8 12], 'XTick', [0:20:120], 'TickDir','out', 'LineWidth', 1.5,'FontName', 'Arial')
end
function dataA = dataIn(MO, surf, force)
%read in mature and pre
if nargin < 1
force = 0;
end
if nargin < 2
showplot = 0;
end
if ispc
indir = fullfile([getenv('HOMEDRIVE') getenv('HOMEPATH') ], 'Dropbox', 'NPCMaturation',...
'expdata', 'emdata');
......@@ -42,12 +309,11 @@ classdef absDelayedMaturation < handle
end
MO.indir = indir;
MO.outdir = outdir;
matfile = fullfile(indir,'emdata.mat');
matfile = fullfile(MO.indir,['emdata' surf '.mat']);
if exist(matfile) && ~force
load(matfile)
return
end
% first entry is minipore second entry is mature pore
data = struct('time', [], 'd_noncore', [], 'd_core_i', [], 'd_core_o', [], 't_core_i', [], 't_core_o', []);
dataA = [data; data];
......@@ -57,7 +323,7 @@ classdef absDelayedMaturation < handle
dataA(2).t_core_i = xlsread(xlsfile, 'Onlyvalues', 'F4:G15');
dataA(2).t_core_o = xlsread(xlsfile, 'Onlyvalues', 'H4:I15');
dataA(2).d_noncore = xlsread(xlsfile, 'Onlyvalues', 'D20:E31');
dataA(2).d_core_i = xlsread(xlsfile, 'Onlyvalues', 'F20:F31');
dataA(2).d_core_i = xlsread(xlsfile, 'Onlyvalues', 'F20:G31');
dataA(2).d_core_o = xlsread(xlsfile, 'Onlyvalues', 'H20:I31');
dataA(1).time = xlsread(xlsfile, 'Onlyvalues', 'K4:K15');
......@@ -69,15 +335,46 @@ classdef absDelayedMaturation < handle
dataA(1).d_core_o = xlsread(xlsfile, 'Onlyvalues', 'P20:Q31');
%% read and create surface area
surfdatafile = fullfile(indir, 'surfaceNucleus_24h_new.txt');
surfdata = load(surfdatafile);
dataA(1).surf = [surfdata(:,1)*60 surfdata(:,2:3)];
dataA(1).surf = [0 surfdata(1,2:3); dataA(1).surf];
% perform a spline fit to smooth the data
dataA(1).surfspl = surfSpline( dataA(1).surf(:,1), dataA(1).surf(:,2),7, 3, 8);
dataA(1).dsurfspl = fnder(dataA(1).surfspl, 1);
switch surf
case '24h'
surfdatafile = fullfile(MO.indir, 'surfaceNucleus_24h_new.txt');
surfdata = load(surfdatafile);
dataA(1).surf = [surfdata(:,1)*60 surfdata(:,2:3)];
dataA(1).surf = [0 surfdata(1,2:3); dataA(1).surf];
% perform a spline fit to smooth the data
dataA(1).surfspl = surfSpline( dataA(1).surf(:,1), dataA(1).surf(:,2),7, 3, 8);
dataA(1).dsurfspl = fnder(dataA(1).surfspl, 1);
dataA(1).surfPar = load(fullfile(MO.indir, 'surfaceNucleusFit_24h_new.txt'));
case '2h'
surfdatafile = fullfile(MO.indir, 'surfaceNucleus_2h.txt');
surfdata = load(surfdatafile);
dataA(1).surf = [surfdata(2:end,1) surfdata(2:end,2:3)];
dataA(1).surf = [[0:1:surfdata(2,1)-1]' repmat(surfdata(2,2:3),7,1) ; dataA(1).surf];
% perform a spline fit to smooth the data
dataA(1).surfspl = surfSpline(dataA(1).surf(:,1), dataA(1).surf(:,2),15, 3, 8);
dataA(1).dsurfspl = fnder(dataA(1).surfspl, 1);
dataA(1).surfPar = load(fullfile(MO.indir, 'surfaceNucleusFit_2h.txt'));
otherwise
error('Usage class(surf, force, showplot), surf= 2h or 24h string')
end
save(matfile, 'dataA');
end
function [A, dAdt, kg] = getSurf(MO, time)
if MO.useSpline
A = fnval(MO.dataF(1).surfspl, time);
dAdt = fnval(MO.dataF(1).dsurfspl, time);
kg = dAdt./A;
else
par = struct('a0', MO.dataA(1).surfPar(1), 'a1', MO.dataA(1).surfPar(2), ...
'kg1', MO.dataA(1).surfPar(3), 'kg2', MO.dataA(1).surfPar(4), 'ts', MO.dataA(1).surfPar(5));
A = par.a0 + par.a1*(1-exp(-par.kg1*(time-par.ts))) + par.kg2*(time-par.ts);
dAdt = (par.a1*par.kg1*exp(-par.kg1*time)+par.kg2).*(time>=par.ts);
kg = dAdt./A;
end
end
end
end
matlabcode/delayedEquationModel/afitround.m 0 → 100644
View file @ 067d60d4
function [dist, distI, distM] = afitround(MO, par, ifit, nrfits, parfile)
denfile = fullfile(MO.outdir,[parfile '_den.txt']);
prodratefile = fullfile(MO.outdir, [parfile '_prodrate.txt']);
fractionMaturefile = fullfile(MO.outdir, [parfile '_fracM.txt']);
parfile = [parfile '.txt'];
MO.runFit(par, ifit, nrfits, 0.5, parfile)
MO.plotBestfit(fullfile(MO.outdir, parfile), 1, 1);
parf = MO.getBestfit(fullfile(MO.outdir, parfile));
[dist, distI, distM] = MO.distMO(parf(ifit), ifit, parf);
dist = sum(dist.^2);
distI = sum(distI.^2);
distM = [sum(distM.^2) sum(distM(1:12).^2) sum(distM(13:24).^2)];
MO.writeDensities(parf,denfile);
MO.writeProdRate(parf, prodratefile);
MO.writeFractionMatureIP(parf, fractionMaturefile);
%%
[soli, solo] = MO.solve(parf, [parf(1) 1200]);
timPM = [180:1:1200];
yi = deval(soli, timPM);
yo = deval(solo, timPM);
A = MO.getSurf(timPM);
mean(yi(3,:)./A*0.68 + yo(3,:)./A*0.32)
std(yi(3,:)./A*0.68 + yo(3,:)./A*0.32)
mean((yi(1,:)+yi(2,:))./A*0.68 + (yi(1,:)+yi(2,:))./A*0.32)
std((yi(1,:)+yi(2,:))./A*0.68 + (yi(1,:)+yi(2,:))./A*0.32)
end
\ No newline at end of file
matlabcode/delayedEquationModel/delayedMaturation.m 0 → 100644
View file @ 067d60d4
%%
% dealyedMaturation model with production proportional to free area
% V = v*A
% par(1): ts, start of interphase assembly
% par(2): v, maximal rate of interphase assembly
% par(3): kM, initiation of maturation
% par(4): tauM, maturation time
% par(5): kd, decay of Mature pore
% par(6): number of minipores at ts inner core
% par(7): number of mature pores at ts inner core
% par(7): number of minipores at ts outer core
% par(8): number of mature pores at ts outer core
classdef delayedMaturation < absDelayedMaturation
methods
function MO = delayedMaturation(surf, useSpline, force)
MO@absDelayedMaturation(surf, useSpline, force)
MO.dataF = MO.dataA;
for i =1:2
MO.dataF(i).noncore = MO.dataF(i).d_noncore;
MO.dataF(i).core_i = MO.dataF(i).d_core_i;
MO.dataF(i).core_o = MO.dataF(i).d_core_o;
end
MO.parnames = 'ts v kM tauM kd P01 M01 P02 M02';
MO.LB = [0, 0, 0.01, 10, 0, 0, 0,0,0]; %Low boundary values
MO.HB = [15, 1, 10, 70, 1, 10000, 10000, 10000, 10000]; %High boundary values
MO.name = mfilename('class')
MO.N = 3;
end
%%
% solve equations for one set of parameters and provide the
function [soli, solo] = solve(MO, par, timeint)
soli = dde23(@MO.dde,[par(4)],@MO.ddehist1, timeint, ddeset, par);
solo = dde23(@MO.dde,[par(4)],@MO.ddehist2, timeint, ddeset, par);
end
%%
% get production rate of minipores and mature pores coming from PM
% only
function [prodRate, maturePM] = getProdRate(MO, par, timeint)
A = MO.getSurf(timeint);
prodRate = [par(2)*(timeint>=par(1)) par(2)*(timeint>=par(1))];
maturePM = [par(7)*exp(-par(5)*(timeint - par(1)))./A par(9)*exp(-par(5)*(timeint - par(1)))./A];
end
function obs = getObservable(MO,sol)
obs(:,1) = sol(1,:)'+ sol(2,:)';
obs(:,2) = sol(3,:)';
end
function s = ddehist1(MO,t,par)
% Constant history function for DDEX1.
s = [0, 0, par(7)];
if t == par(1)
s = [par(6), 0, par(7)];
end
end
function s = ddehist2(MO,t,par)
% Constant history function for DDEX1.
s = [0, 0, par(9)];
if t == par(1)
s = [par(8), 0, par(9)];
end
end
function dydt = dde(MO, t, y, Z, par)
% par(1): ts, start of interphase assembly
% par(2): v, maximal rate of interphase assembly
% par(3): kM, initiation of maturation
% par(4): tauM, maturation time
% par(5): kd, decay of Mature pore
% par(6): number of minipores at ts
% par(7): number of mature pores at 0
% Delay Differential equations function for DDEX1.
ylag = Z(:,1);
if t > par(1)
v = par(2);
else
v = 0;
end
A = MO.getSurf(t);
dydt = [ v*A - par(3)*y(1)
par(3)*y(1) - par(3)*ylag(1)
par(3)*ylag(1) - par(5)*y(3)];
end
function [par, norm] = testFit(MO)
ts = 10;
At0 = MO.getSurf(ts);
P01 = 4.493*At0;
M01 = 4.79*At0;
P02 = 3.935*At0;
M02 = 9.048*At0;
if nargin < 2
par = [ts, 0.014, 1.136, 40, 4.2e-04, P01, M01, P02, M02];
end
ifit = [2:4 6:9];
[parf, norm] = MO.fitModel(ifit, par)
par(ifit) = parf;
MO.plotResult(par)
end
end
end
\ No newline at end of file
matlabcode/delayedEquationModel/delayedMaturationDistributedDelay.m 0 → 100644
View file @ 067d60d4
%%
% dealyedMaturation model with production proportional to free area
% V = v*A.
% delay is uniform distributed along tauM1 and tauM2
% par(1): ts, start of interphase assembly
% par(2): v, maximal rate of interphase assembly
% par(3): kM, initiation of maturation
% par(4): tauM1, lowest maturation time
% par(5): DtauM, tauM1 + DtauM gives the highest maturation time tauM2
% par(6): kd, decay of Mature pore
% par(7): number of minipores at ts inner core
% par(8): number of mature pores at ts inner core
% par(9): number of minipores at ts outer core
% par(10): number of mature pores at ts outer core
classdef delayedMaturationDistributedDelay < absDelayedMaturation
methods
function MO = delayedMaturationDistributedDelay(surf, useSpline, force)
MO@absDelayedMaturation(surf, useSpline, force)
MO.dataF = MO.dataA;
for i =1:2
MO.dataF(i).noncore = MO.dataF(i).d_noncore;
MO.dataF(i).core_i = MO.dataF(i).d_core_i;
MO.dataF(i).core_o = MO.dataF(i).d_core_o;