Skip to content

GitLab

Commit 53505883 authored by Antonio Politi's avatar Antonio Politi
Browse files

first commit of expdata and some code

parents
Expand all Hide whitespace changes
Inline Side-by-side
expdata/livecell/README.txt 0 → 100644
View file @ 53505883
HeLa4D.xls contains the analyzed innercore and non-core intensities for dividing cells
2xZFN mEGFP-Nup107 26, 31
and
CRISPR mEGFP-Nup358/RanBP2 118
from file HeLaNup358-2h-background-individually-average-intensity-forAntonio.xls
from file HeLaNup107-2h-160121-background-individual-average-intensity-forAntonio.xls
data acquired by Shotaro Otsuka
expdata/livecell_siRNA/README.txt 0 → 100644
View file @ 53505883
Pom121_siRNA.xls contains the analyzed innercore and non-core intensities for dividing cells
CRISPR mEGFP-Nup358/RanBP2 118
data acquired by Joe Padget
Pom121 has been knocked-down for 48h (?)
matlabcode/delayedEquationModel/delayedMaturation.m 0 → 100644
View file @ 53505883
This diff is collapsed.
matlabcode/odeModel/absPoreMaturation.m 0 → 100644
View file @ 53505883
%%
% Abstact class for model of pore maturation
%%
classdef absPoreMaturation < handle
properties (Access = protected)
timeI = [0:150]; %intergration time points
protA; %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]
end
properties
N1; %number of steps controlled by rate constant k1
N2; %number of steps controlled by rate constant k2
indir;
outdir;
protF; %data structure containing data to be fitted
end
methods (Abstract)
A = modelMatrix(MO,par);
dydt = model(MO, par);
prot = computeProteins(MO, par, iprot)
end
methods
function MO = absPoreMaturation(N1, N2, tpm, ti, force, showplot)
if nargin < 5
force = 0;
end
if nargin < 6
showplot = 0;
end
MO.setTimeStart(tpm, ti);
MO.N1 = N1;
MO.N2 = N2;
MO.protA = MO.dataIn(force, showplot);
end
function setTimeStart(MO, tpm, ti)
% set time of start of assembly of postmitotic and interphase
assert(size(tpm) > 1)
assert(size(ti) > 1)
MO.tpm = round(tpm);
MO.ti = round(ti);
end
function [idxm, idxd] = findIdxs(MO)
idxm = [1 1];
idxd = [1 1];
% for i=1:2
% it_tmp = find(MO.protF(1).time == MO.tpm(i));
% idxm(i) = 1;
% if isempty(it_tmp)
% idxm(i) = find(MO.timeI(MO.tpm(i)+1:end) == MO.protF(1).time(1));
% idxd(i) = 1;
% else
% idxd(i) = it_tmp;
% end
% end
end
function parloc = getparset(MO, parfit, ifit, par)
parloc = par;
parloc(ifit) = parfit;
end
function dist = distData(MO, parfit, ifit, par)
%% calculate distance model to data
parloc = MO.getparset(parfit, ifit, par);
prot = MO.computeProteins(parloc, [1 2]);
% if useStd == 0 do not divide by standard deviation
[idxm, idxd] = MO.findIdxs();
useStd = 1;
distNC = [];
distC = [];
for i=1:2
if useStd
distNC = [distNC (MO.protF(i).noncore(idxd(1):end,1)' - prot(i).NC(idxm(1):end))./MO.protF(i).noncore(idxd(1):end,2)'];
distC = [distC (MO.protF(i).core(idxd(2):end,1)' - prot(i).C(idxm(2):end))./MO.protF(i).core(idxd(2):end,2)'];
else
distNC = [distNC MO.protF(i).noncore(idxd(1):end,1)' - prot(i).NC(idxm(1):end)];
distC = [distC MO.protF(i).core(idxd(2):end,1)' - prot(i).C(idxm(2):end)];
end
end
dist = [distNC distC];
end
%%
% function to create a graph of the data
function showGraph(MO, par, iprot)
prot = MO.computeProteins(par, iprot);
[idxm, idxd] = MO.findIdxs()
for i=1:2
setupFigure(i, [200 200 500 700]);
clf
subplot(3,1,1)
title(MO.protF(i).name)
hold
plot(MO.protF(i).time, MO.protF(i).noncore(:,1), 'bo');
plot(MO.protF(i).time, MO.protF(i).core(:,1), 'go');
plot(MO.protF(i).time(idxd(1):end), prot(i).NC(idxm(1):end), 'b--');
plot(MO.protF(i).time(idxd(2):end), prot(i).C(idxm(2):end), 'g--');
ylabel('Relative fluorescence')
legend('Core data', 'Non-core data', 'Core model', 'Non-core model')
ylim([0 1.2])
subplot(3,1,2)
hold
plot(MO.protF(i).time, MO.protF(i).noncore(:,1), 'bo');
plot(MO.protF(i).time(idxd(1):end), prot(i).NC_pm(idxm(1):end)', 'b--');
plot(MO.protF(1).time(idxd(1):end), prot(i).NC_ip(idxm(1):end)', 'b--', 'LineWidth', 2);
ylabel('Relative fluorescence')
legend('Non-Core data', 'Non-Core 2', 'Non-Core 2')
ylim([0 1.2])
subplot(3,1,3)
hold
plot(MO.protF(i).time, MO.protF(i).core(:,1), 'go');
plot(MO.protF(i).time(idxd(2):end), prot(i).C_pm(idxm(2):end)', 'g--');
plot(MO.protF(i).time(idxd(2):end), prot(i).C_ip(idxm(2):end)', 'g--', 'LineWidth', 2);
legend('Core data', 'Core 1', 'Core 2')
ylabel('Relative fluorescence')
xlabel('Time (min)')
ylim([0 1.2])
end
end
function y = solveModelN(MO, par, time, iprot)
y = ode45(@MO.model,time,[1;zeros(MO.N1+MO.N2,1)], odeset, par, iprot);
end
%%
% read in the experimental data
%
function protA = dataIn(MO, force, showplot)
if nargin < 1
force = 0;
end
if nargin < 2
showplot = 1;
end
if ispc
indir = fullfile([getenv('HOMEDRIVE') getenv('HOMEPATH') ], 'Dropbox', 'Shotaro_PoreMaturation',...
'ExpData', 'LiveCellImaging');
outdir = fullfile([getenv('HOMEDRIVE') getenv('HOMEPATH') ], 'Dropbox', 'Shotaro_PoreMaturation',...
'Simulations', 'LiveCellImaging');
end
if ismac
indir = fullfile(getenv('HOME'), 'Dropbox', 'Shotaro_PoreMaturation',...
'ExpData', 'LiveCellImaging');
outdir = fullfile(getenv('HOME'), 'Dropbox', 'Shotaro_PoreMaturation',...
'Simulations', 'LiveCellImaging');
end
MO.indir = indir;
MO.outdir = outdir;
matfile = fullfile(indir,'livecelldata.mat');
%load matfile if matfile exists or reading is not forced
if exist(matfile) && ~force
load(matfile)
MO.timeI = [0:protA(1).time(end)]
return
end
% struct containing the experimental data
prot = struct('name','Nup107', 'time',[], 'd_core',[], 'd_noncore', [], ...
'd_core_cells', [],'d_noncore_cells', [], 'tot_core', [], 'tot_noncore', [])
protA = [prot;prot]
protA(2).name = 'Nup358';
% location of XLS files
xlsfile = fullfile(indir, 'HeLa4D.xls');
sheet = {'Nup107_AVG', 'Nup358_AVG'};
%% Raed Nup107 data
protA(1).d_core_cells = xlsread(xlsfile, 'Nup107_AVG', 'C11:AF132');
protA(1).d_noncore_cells = xlsread(xlsfile, 'Nup107_AVG', 'AH11:BK132');
protA(1).d_core = [xlsread(xlsfile, 'Nup107_AVG', 'BN11:BN132') ...
xlsread(xlsfile, 'Nup107_AVG', 'BP11:BP132')];
protA(1).d_noncore = [xlsread(xlsfile, 'Nup107_AVG', 'BO11:BO132') ...
xlsread(xlsfile, 'Nup107_AVG', 'BQ11:BQ132')];
protA(1).time = xlsread(xlsfile, 'Nup107_AVG', 'B11:B132');
protA(1).tot_core = [xlsread(xlsfile, 'Nup107_total', 'BN11:BN132') xlsread(xlsfile, 'Nup107_total', 'BP11:BP132')];
protA(1).tot_noncore = [xlsread(xlsfile, 'Nup107_total', 'BO11:BO132') xlsread(xlsfile, 'Nup107_total', 'BQ11:BQ132')];
protA(2).d_core_cells = xlsread(xlsfile, 'Nup358_AVG', 'C11:AA132');
protA(2).d_noncore_cells = xlsread(xlsfile, 'Nup358_AVG', 'AC11:BA132');
protA(2).d_core = [xlsread(xlsfile, 'Nup358_AVG', 'BD11:BD132') ...
xlsread(xlsfile, 'Nup358_AVG', 'BF11:BF132')];
protA(2).d_noncore = [xlsread(xlsfile, 'Nup358_AVG', 'BE11:BE132') ...
xlsread(xlsfile, 'Nup358_AVG', 'BG11:BG132')];
protA(2).time = xlsread(xlsfile, 'Nup358_total', 'B11:B132');
protA(2).tot_core = [xlsread(xlsfile, 'Nup358_total', 'BD11:BD132') xlsread(xlsfile, 'Nup358_total', 'BF11:BF132')];
protA(2).tot_noncore = [xlsread(xlsfile, 'Nup358_total', 'BE11:BE132') xlsread(xlsfile, 'Nup358_total', 'BG11:BG132')];
%% read surface area
protA(1).surf = xlsread(xlsfile, 'NUP107_area', 'AH11:AI132');
protA(2).surf = xlsread(xlsfile, 'NUP358_area', 'AC11:AD132');
% perform a spline fit to smooth the data
protA(1).surfspl = surfSpline(protA(1).time, protA(1).surf(:,1),7, 3, 8);
protA(2).surfspl = surfSpline(protA(2).time, protA(2).surf(:,1),7, 3, 8);
% compute derivative of spline fit to be used when solving the
% model that includes the area
protA(1).dsurfspl = fnder(protA(1).surfspl, 1);
protA(2).dsurfspl = fnder(protA(2).surfspl, 1);
protA(1).tot_core2 = protA(1).d_core.*repmat(fnval(protA(1).surfspl, protA(1).time),1,2);
protA(1).tot_noncore2 = protA(1).d_noncore.*repmat(fnval(protA(1).surfspl, protA(1).time),1,2);
protA(2).tot_core2 = protA(2).d_core.*repmat(fnval(protA(2).surfspl, protA(2).time),1,2);
protA(2).tot_noncore2 = protA(2).d_noncore.*repmat(fnval(protA(2).surfspl, protA(2).time),1,2);
for i=1:2
sP = sort(protA(i).tot_core2(:,1));
protA(i).tot_core2 = protA(i).tot_core2/mean(sP(end-3:end));
sP = sort(protA(i).tot_noncore2(:,1));
protA(i).tot_noncore2 = protA(i).tot_noncore2/mean(sP(end-3:end));
end
save(matfile, 'protA');
%%
if showplot
figure(1)
clf
col = {'b', 'r'};
for i=1:length(protA)
subplot(1,2,i)
hold
errorbar(protA(i).time, protA(i).tot_core(:,1), protA(i).d_core(:,2), 'g')
errorbar(protA(i).time, protA(i).tot_noncore(:,1), protA(i).tot_noncore(:,2), 'b')
end
end
end
end
end
matlabcode/odeModel/findCI_poreMaturation.m 0 → 100644
View file @ 53505883
%%
% find confidence Interval using likelihood profile method for ode model.
% The scaling coefficients are not varied
% pM: class of model
% filenames: 4 filenames containing parameter values
function findCI_poreMaturation(pM, filenames)
try
k1_pm = load(filenames{1});
k2_pm = load(filenames{2});
k1_ip = load(filenames{3});
k2_ip = load(filenames{4});
catch
%priming
tpm = [4 4];
% time begin interphase assembly for Non-core and core region
tip = [10 10];
hb = [2 2 2 2 2 2 2 2];
lb = [0 0 0 0 0 0 0 0];
opti = optimset('Display', 'iter');
solB = [0.3549 0.0449 0.0435 0.0492 0.92 0.5 1 1 1 1]; %bset solution
ifit = [1 2 3 4]; %parameters to fit
for i=1:4
iv = hb.*rand(1,length(hb))
[parFit{i}, norm(i)] = lsqnonlin(@pM.distData,iv, lb, hb, opti, [1:4], solB);
end
normB = min(norm);
ibf = find(norm == normB);
solB = pM.getparset(parFit{ibf}, ifit, solB);
%% vary each one parameter
minP = 0.2;
fitParCI = {[2 3 4], [1 3 4], [1 2 4], [1 2 3]};
%fit parameter [2 3 4] and vary parameter 1 step by step, etc.
step_p = 100;
for i=1:4
hb = [2 2 2];
lb = [0 0 0];
fidloc = fopen(filenames{i}, 'w');
fprintf(fidloc, '%%k1_pm \t k2_pm \t k1_ip \t k2_ip \t norm \t normBest\n');
ifit = fitParCI{i};
solLoc = solB;
parVar = solB(i)*10.^((-minP+minP*2*[1:step_p]/step_p));
for j = 1:step_p
solLoc(i) = parVar(j);
iv = solLoc(ifit);
[parFit, norm] = lsqnonlin(@pM.distData,iv, lb, hb, opti, ifit, solLoc);
sol = pM.getparset(parFit, ifit, solLoc);
fprintf(fidloc, '%.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e\n', sol(1), sol(2), sol(3), sol(4), norm, normB);
end
fclose(fidloc);
end
k1_pm = load(filenames{1});
k2_pm = load(filenames{2});
k1_ip = load(filenames{3});
k2_ip = load(filenames{4});
end
%%
NrDataPts = size(pM.protF(1).core,1)+size(pM.protF(2).core,1) + size(pM.protF(1).noncore,1)+size(pM.protF(2).noncore,1);
parList = {k1_pm, k2_pm, k1_ip, k2_ip};
for i=1:4
par = parList{i};
CI = par(:,i).*((abs(NrDataPts*log(par(:,end-1)/NrDataPts) - NrDataPts*log(par(:,end)/NrDataPts)))<chi2inv(0.95,1));
CIL(i,:) = [par(find(par(:,end-1)-par(:,end)<1e-5), i) min(CI(find(CI>0))) max(CI(find(CI>0)))];
end
end
matlabcode/odeModel/findCI_poreMaturation_coeffMax.m 0 → 100644
View file @ 53505883
%%
% find confidence Interval using likelihood profile method for ode model.
% The scaling coefficients are varied
% pM: Class of model
% filenames: 4 filenames containing parameter values
function findCI_poreMaturation_coeffMax(pM, filenames)
try
k1_pm = load(filenames{1});
k2_pm = load(filenames{2});
k1_ip = load(filenames{3});
k2_ip = load(filenames{4});
catch
%priming
tpm = [4 4];
% time begin interphase assembly for Non-core and core region
tip = [10 10];
hb = [2 2 2 2 2 2 2 2];
lb = [0 0 0 0 0 0 0 0];
opti = optimset('Display', 'iter');
solB = [0.3549 0.0449 0.0435 0.0492 0.92 0.5 1 1 1 1]; %bset solution
ifit = [1 2 3 4 7:10]; %parameters to fit
for i=1:4
iv = hb.*rand(1,length(hb))
[parFit{i}, norm(i)] = lsqnonlin(@pM.distData,iv, lb, hb, opti, ifit, solB);
end
normB = min(norm);
ibf = find(norm == normB);
solB = pM.getparset(parFit{ibf}, ifit, solB);
%% vary each one parameter
minP = 0.2;
fitParCI = {[2 3 4 7:10], [1 3 4 7:10], [1 2 4 7:10], [1 2 3 7:10]};
step_p = 100;
for i=1:4
hb = [2 2 2 1.2 1.2 1.2 1.2];
lb = [0 0 0 0 0 0 0];
fidloc = fopen(filenames{i}, 'w');
fprintf(fidloc, '%%k1_pm \t k2_pm \t k1_ip \t k2_ip \t norm \t normBest\n');
ifit = fitParCI{i};
solLoc = solB;
parVar = solB(i)*10.^((-minP+minP*2*[1:step_p]/step_p));
for j = 1:step_p
solLoc(i) = parVar(j);
iv = solLoc(ifit);
[parFit, norm] = lsqnonlin(@pM.distData,iv, lb, hb, opti, ifit, solLoc);
sol = pM.getparset(parFit, ifit, solLoc);
fprintf(fidloc, '%.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \n', sol(1), sol(2), sol(3), ...
sol(4), sol(7), sol(8), sol(9), sol(10), norm, normB);
end
fclose(fidloc);
end
k1_pm = load(filenames{1});
k2_pm = load(filenames{2});
k1_ip = load(filenames{3});
k2_ip = load(filenames{4});
end
%%
NrDataPts = size(pM.protF(1).core,1)+size(pM.protF(2).core,1) + size(pM.protF(1).noncore,1)+size(pM.protF(2).noncore,1);
parList = {k1_pm, k2_pm, k1_ip, k2_ip};
for i=1:4
par = parList{i};
CI = par(:,i).*((abs(NrDataPts*log(par(:,end-1)/NrDataPts) - NrDataPts*log(par(:,end)/NrDataPts)))<chi2inv(0.95,1));
CIL(i,:) = [par(50, i) min(CI(find(CI>0))) max(CI(find(CI>0)))];
end
CIL
end
% %systematically change tauM and kM with fixed parameters (do not refit)
% function findCI(par, fileout1, fileout2)
% try
% kM = load(fileout1);
% tauM = load(fileout2);
% catch
% stepi = 50;
% stepj = 50;
% parloc = par;
% data = getData(density, 0, 0);
% %options = odeset('MaxStep', 0.1);
% clear('normFit')
% minTMv = 0.2;
% minkMV = 1.6;
%
% tauMv = par.tauM*10.^((-minTMv+minTMv*2*[1:stepi]/stepi));
% kMv = par.kM*10.^((-minkMV + minkMV*2*[1:stepj]/stepj));
% tofit_names_loc = {'tauM', 'P01', 'M01', 'P02', 'M02'};
% normO = sum(distModel(getfields(par, tofit_names_loc), tofit_names_loc, par, data, options).^2);
%
% tofit_names_loc = {'tauM', 'P01', 'M01', 'P02', 'M02'};
% fidloc = fopen(fileout1, 'w');
% fprintf(fidloc, '%%kM \t tauM \t P01 \t M01 \t P02 \t M02 \t norm \t normBest\n');
% parloc = par;
% for idx = 1:stepj
% parloc.tauM = par.tauM;
% parloc.kM = kMv(idx);
% [parfit_kM, norm_kM(idx)] = performFit(tofit_names_loc, lbStruct, hbStruct, parloc, 5, 0.1)
% parloc = setfields(parloc, tofit_names_loc, parfit_kM);
% fprintf(fidloc, '%.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e\n', parloc.kM, parloc.tauM, parloc.P01, parloc.M01, parloc.P02, parloc.M02, norm_kM(idx), normO);
% end
% fclose(fidloc);
%
% tofit_names_loc = {'kM', 'P01', 'M01', 'P02', 'M02'};
% fidloc = fopen(fileout2, 'w');
% fprintf(fidloc, '%%kM \t tauM \t P01 \t M01 \t P02 \t M02 \t norm \t normBest\n');
% parloc = par;
% for idx = 1:31
% parloc.tauM = tauMv(idx);
% parloc.kM = 2;
% [parfit_tauM, norm_tauM(idx)] = performFit(tofit_names_loc, lbStruct, hbStruct, parloc, 5, 0.1)
% parloc = setfields(parloc, tofit_names_loc, parfit_tauM);
% fprintf(fidloc, '%.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e \t %.3e\n', parloc.kM, parloc.tauM, parloc.P01, parloc.M01, parloc.P02, parloc.M02, norm_tauM(idx), normO);
% end
% fclose(fidloc);
% end
% %%
% data = getData(density, 0, 0);
% NrDataPts = (size(data,1)-1)/2*size(data,2);
% CIkM = kM(:,1).*((abs(NrDataPts*log(kM(:,end-1)/NrDataPts) - NrDataPts*log(kM(:,end)/NrDataPts)))<chi2inv(0.95,1));
% CIkM = [min(CIkM(find(CIkM>0))) max(CIkM(find(CIkM>0)))]
% CItauM = tauM(:,2).*((abs(NrDataPts*log(tauM(:,end-1)/NrDataPts) - NrDataPts*log(tauM(:,end)/NrDataPts)))<chi2inv(0.95,1));
% CItauM = [min(CItauM(find(CItauM>0))) max(CItauM(find(CItauM>0)))]
% end
\ No newline at end of file
matlabcode/odeModel/poreMaturationDen.m 0 → 100644
View file @ 53505883
classdef poreMaturationDen < absPoreMaturation
%% Multi step maturation process N1 steps with rate equal k1 and N2 steps with rate equal k2
% Use surface area to compute density
% k1 set the rate for formation of protein1 containing complex
% k2 set the rate of intermediate steps until protein2 is
% integrated
% X_1 -k1> X_2 ... -k1> X_(N1+1) -k2> X_(N1+2) ... -k2> X_(N1+N2+1)
% sum X_(i=N1+1)_end: protein1
% X_(N1+N2+1): protein2
methods
function MO = poreMaturationDen(N1, N2, tpm, ti, force, showplot)
MO@absPoreMaturation(N1, N2, tpm, ti,force, showplot)
MO.protF = MO.protA;
for i=1:2
MO.protF(i).core = MO.protA(i).d_core;
MO.protF(i).noncore = MO.protA(i).d_noncore;
end
end
function setTimeStart(MO, tpm, ti)
% set time of start of assembly of postmitotic and interphase
assert(length(tpm) > 1)
assert(length(ti) > 1)
MO.tpm = round(tpm);
MO.ti = round(ti);
end
function [prot] = computeProteins(MO, par,iprot, varargin)
%% compute distance model to data
% par: par(1) - k1 for postmitotic
% par(2) - k2 for postmitotic
% par(3) - k1 for interphase
% par(4) - k2 for interphase
% par(5) - fraction postmitotic in non-core region
% par(6) - fraction postmitotic in core region
% par(7) - coefficient scaling non-core protein 1
% par(8) - coefficient scaling non-core protein 2
% par(9) - coefficient scaling core protein 1
% par(10) - coefficient scaling core protein 2
if nargin <= 3
tend = 125;
else
tend = varargin{1};
end
%experimental data has 1 min time resolution
p1_NC_pm = MO.solveModelN([par(1) par(2)], [MO.tpm(1) tend],iprot(1));
p1_NC_ip = MO.solveModelN([par(3) par(4)], [MO.ti(1) tend],iprot(1));
p2_NC_pm = MO.solveModelN([par(1) par(2)], [MO.tpm(1) tend],iprot(2));
p2_NC_ip = MO.solveModelN([par(3) par(4)], [MO.ti(1) tend],iprot(2));
p1_C_pm = MO.solveModelN([par(1) par(2)], [MO.tpm(2) tend],iprot(1));
p1_C_ip = MO.solveModelN([par(3) par(4)], [MO.ti(2) tend],iprot(1));
p2_C_pm = MO.solveModelN([par(1) par(2)], [MO.tpm(2) tend],iprot(2));
p2_C_ip = MO.solveModelN([par(3) par(4)], [MO.ti(2) tend],iprot(2));
pm = deval(p1_NC_pm,[MO.tpm(1):1:tend]);
ip = deval(p1_NC_ip,[MO.ti(1):1:tend]);
p1NC_pm = [zeros(1,MO.tpm(1)-4) sum(pm(MO.N1+1:end,:))]*par(5)*par(7);
p1NC_ip = [zeros(1,MO.ti(1)-4) sum(ip(MO.N1+1:end,:))]*(1-par(5))*par(7);
p1NC = (p1NC_pm+p1NC_ip);
pm = deval(p2_NC_pm,[MO.tpm(1):1:tend]);
ip = deval(p2_NC_ip,[MO.ti(1):1:tend]);
p2NC_pm = [zeros(1,MO.tpm(1)-4) pm(end,:)]*par(5)*par(8);
p2NC_ip = [zeros(1,MO.ti(1)-4) ip(end,:)]*(1-par(5))*par(8);
p2NC = (p2NC_pm + p2NC_ip);
pm = deval(p1_C_pm,[MO.tpm(2):1:tend]);
ip = deval(p1_C_ip,[MO.ti(2):1:tend]);
p1C_pm = [zeros(1,MO.tpm(2)-4) sum(pm(MO.N1+1:end,:))]*par(6)*par(9);
p1C_ip = [zeros(1,MO.ti(2)-4) sum(ip(MO.N1+1:end,:))]*(1-par(6))*par(9);
p1C = (p1C_pm + p1C_ip);
pm = deval(p2_C_pm,[MO.tpm(2):1:tend]);
ip = deval(p2_C_ip,[MO.ti(2):1:tend]);
p2C_pm = [zeros(1,MO.tpm(2)-4) pm(end,:)]*par(6)*par(10);
p2C_ip = [zeros(1,MO.ti(2)-4) ip(end,:)]*(1-par(6))*par(10);
p2C = (p2C_pm + p2C_ip);
p1 = struct('time_NC', [4:tend],'time_C', [4:tend], 'NC', p1NC, 'NC_pm', p1NC_pm, 'NC_ip', p1NC_ip,...
'C', p1C, 'C_pm', p1C_pm, 'C_ip', p1C_ip, 'C_pm_pure', p1C_pm/par(6)/par(9), 'C_ip_pure', p1C_ip/(1-par(6))/par(9));
p2 = struct('time_NC',[4:tend],'time_C',[4:tend], 'NC', p2NC, 'NC_pm', p2NC_pm, 'NC_ip', p2NC_ip, ...
'C', p2C, 'C_pm', p2C_pm, 'C_ip', p2C_ip, 'C_pm_pure', p2C_pm/par(6)/par(10), 'C_ip_pure', p2C_ip/(1-par(6))/par(10));
prot = [p1;p2];
end
function [idxm, idxd] = findIdxs(MO)
idxm = [1 1];
idxd = [1 1];
end
function A = modelMatrix(MO, par)
%% Generate matrix for multi step maturation process N1 steps with rate equal k1 and N2 steps with rate equal k2