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

Multiparameter live cell model

parent b97eff37
......@@ -83,6 +83,7 @@ classdef absPoreMaturation < handle
end
end
dist = [distNC distC];
%dist = [distNC;
end
%%
......@@ -99,8 +100,8 @@ classdef absPoreMaturation < handle
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--');
plot(MO.protF(i).time(idxd(1):end), prot(i).NC(idxm(1):end), 'c--');
plot(MO.protF(i).time(idxd(2):end), prot(i).C(idxm(2):end), 'r--');
ylabel('Relative fluorescence')
legend('Core data', 'Non-core data', 'Core model', 'Non-core model')
ylim([0 1.2])
......@@ -117,8 +118,8 @@ classdef absPoreMaturation < handle
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);
plot(MO.protF(i).time(idxd(2):end), prot(i).C_pm(idxm(2):end)', 'r--');
plot(MO.protF(i).time(idxd(2):end), prot(i).C_ip(idxm(2):end)', 'r--', 'LineWidth', 2);
legend('Core data', 'Core 1', 'Core 2')
ylabel('Relative fluorescence')
xlabel('Time (min)')
......@@ -126,7 +127,12 @@ classdef absPoreMaturation < handle
end
end
function [A, dAdt, kg] = getSurf(MO, time)
A = fnval(MO.protF(1).surfspl, time);
dAdt = fnval(MO.protF(1).dsurfspl, time);
kg = dAdt./A;
end
function y = solveModelN(MO, par, time, iprot)
y = ode45(@MO.model,time,[1;zeros(MO.N1+MO.N2,1)], odeset, par, iprot);
end
......
%%
% Abstact class for model of pore maturation
%%
classdef absPoreMaturationMultiPar < 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
N; %total number of steps
N1; % position where 1st protein is present
indir; %input directory
outdir; %output directory
LB; %lower boundary
HB; %Higher boundary
protF; %data structure containing data to be fitted
parnames;
name;
end
methods (Abstract)
A = modelMatrix(MO,par);
dydt = model(MO, par);
prot = computeProteins(MO, par, iprot)
end
methods
function MO = absPoreMaturationMultiPar(N, N1, tpm, ti, force, showplot)
if nargin < 5
force = 0;
end
if nargin < 6
showplot = 0;
end
MO.setTimeStart(tpm, ti);
MO.N = N;
MO.N1 = N1;
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
%%
% 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
%%
% 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);
fprintf(fid,'%% N1 %d\n', MO.N1);
parnamesA = strsplit( MO.parnames, ' ');
fprintf(fid, '%% fittedIds %s\n', num2str(ifit));
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 [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 = distMO(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];
%dist = [distNC;
end
function par = getBestfit(MO, filename)
try
parV = load(filename);
catch
error(sprintf('could not load %s\n', filename));
end
[val, ip] = min(parV(:,end));
par = parV(ip, 1:end-1);
end
%%
% function to create a graph of the data
function showGraph(MO, par, iprot, save, filename)
prot = MO.computeProteins(par, iprot);
[idxm, idxd] = MO.findIdxs()
ifig = 1;
setupFigure(ifig, [200 200 700 700]);
clf
for i=1:2
subplot(3,2,1+(i-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), 'c--');
plot(MO.protF(i).time(idxd(2):end), prot(i).C(idxm(2):end), 'r--');
ylabel('Relative fluorescence')
hl = legend('Core data', 'Non-core data', 'Core model', 'Non-core model')
set(hl,'Location', 'SouthEast')
ylim([0 1.2])
subplot(3,2,3+(i-1))
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)', 'c--');
plot(MO.protF(1).time(idxd(1):end), prot(i).NC_ip(idxm(1):end)', 'c--', 'LineWidth', 2);
ylabel('Relative fluorescence')
hl = legend('Non-Core data', 'Non-Core post-mito', 'Non-Core ip')
set(hl,'Location', 'SouthEast')
ylim([0 1.2])
subplot(3,2,5+(i-1))
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)', 'r--');
plot(MO.protF(i).time(idxd(2):end), prot(i).C_ip(idxm(2):end)', 'r--', 'LineWidth', 2);
hl = legend('Core data', 'Core post-mito', 'Core ip')
set(hl,'Location', 'SouthEast')
ylabel('Relative fluorescence')
xlabel('Time (min)')
ylim([0 1.2])
end
if save
[pathstr, filen, exte] = fileparts(filename);
saveas(ifig, fullfile(pathstr, [filen '.pdf']));
end
end
function Den = writeDensities(MO, par, outfile)
fid(1) = fopen(fullfile(MO.outdir, [outfile '_NUP107.csv']), 'w')
fid(2) = fopen(fullfile(MO.outdir, [outfile '_NUP358.csv']), 'w')
for i = 1:2
fprintf(fid(i), 'time\tNC\tNC_pm\tNC_ip\tC\tC_pm\tC_ip\tC_ip_pure\tC_ip_density\tC_ip_density_nop2\tC_ip_density_nop1p2\tC_density\tC_density_nop2\tC_density_nop1p2\n');
prot = MO.computeProteins(par, [1 2])
for it = 1:length(prot(i).time_C)
fprintf(fid(i), '%d\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n', prot(i).time_NC(it), prot(i).NC(it), prot(i).NC_pm(it), prot(i).NC_ip(it), ...
prot(i).C(it), prot(i).C_pm(it), prot(i).C_ip(it), prot(i).C_ip_pure(it), prot(i).C_ip_den(it), prot(i).C_ip_den_nop2(it),prot(i).C_ip_den_nop1p2(it),...
prot(i).C_den(it), prot(i).C_den_nop2(it), prot(i).C_den_nop1p2(it));
end
fclose(fid(i));
end
end
function [A, dAdt, kg] = getSurf(MO, time)
A = fnval(MO.protF(1).surfspl, time);
dAdt = fnval(MO.protF(1).dsurfspl, time);
kg = dAdt./A;
end
function y = solveModelN(MO, par, time, iprot)
y = ode45(@MO.model,time,[1;zeros(MO.N-1,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', 'NPCMaturation',...
'expdata', 'livecell');
outdir = fullfile([getenv('HOMEDRIVE') getenv('HOMEPATH') ], 'Dropbox', 'NPCMaturation',...
'results', 'livecell');
end
if ismac
indir = fullfile(getenv('HOME'), 'Dropbox', 'NPCMaturation',...
'expdata', 'livecell');
outdir = fullfile(getenv('HOME'), 'Dropbox', 'NPCMaturation',...
'results', 'livecell');
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
%% collections of function calls to run the fits and compute confidence intervals
% time begin interphase assembly for Non-core and core region
clear;
clear all;
fclose all
tip = [10 10];
tpm = [4 4];
nrfits = 10;
for N=3:10
for N1 = 2:N-1
MO = poreMaturationTotalProdDeg(N,N1, tpm, tip, 0, 0);
par = [1*ones(1,2*(N-1)) 0 0 0.92 0.5 1 1 1 1];
ifit = [1:2*(N-1) 2*N+1:2*N+6];
fname = ['poreMatuLiveMultiParFractionN_' num2str(N) 'N1_' num2str(N1) '.txt'];
MO.runFit(par, ifit, nrfits, 0.5, fname);
parf = MO.getBestfit(fullfile(MO.outdir, fname));
MO.showGraph(parf,[1 2], 1, fullfile(MO.outdir, fname) );
ifit = [1:2*(N-1) 2*N+3:2*N+6];
fname = ['poreMatuLiveMultiParFractionN_' num2str(N) 'N1_' num2str(N1) '.txt'];
MO.runFit(par, ifit, nrfits, 0.5, fname)
parf = MO.getBestfit(fullfile(MO.outdir, fname));
MO.showGraph(parf,[1 2], 1, fullfile(MO.outdir, fname) );
end
end
%%
clear;
clear all;
fclose all
tip = [10 10];
tpm = [4 4];
nrfits = 10;
for N=3:10
for N1 = 2:N-1
MO = poreMaturationTotalProdDeg(N,N1, tpm, tip, 0, 0);
par = [1*ones(1,2*(N-1)) 4.200e-04 4.2e-4 0.92 0.5 1 1 1 1];
ifit = [1:2*N 2*N+1:2*N+6];
MO.runFit(par, ifit, nrfits, 0.5, ['poreMatuLiveMultiParFractionProdDegN_' num2str(N) 'N1_' num2str(N1) '.txt'])
ifit = [1:2*N 2*N+3:2*N+6];
MO.runFit(par, ifit, nrfits, 0.5, ['poreMatuLiveMultiParProdDegN_' num2str(N) 'N1_' num2str(N1) '.txt'])
end
end
\ No newline at end of file
classdef poreMaturationReducedPara < absPoreMaturationMultiPar
%% Multi step maturation process N1 steps with rate equal k1 and N2 steps with rate equal k2
% Compute total number of pores so no need to use the surface area
% k1-kN-1 are the transition rates for postmitotic and interphase. The
% last transition is the accumulation of the 2nd protein (NUP358)
% N1 is the position where the first protein is present
% X_1 -k1> X_2 ... -k2> X_(N1+1) -k3> X_(N1+2) ... -kN-1> X_(N)
% sum X_(i=N1)_end: protein1
% X_(N): protein2
% paramters are
% k1, k2 (postmitotic); k1, k2 (interphase); kd; vproduction;
% fraction_maturePore_noncore; fraction_maturePore_core;
% scaling_factor_protein1_NC;scaling_factor_protein2_NC;scaling_factor_protein1_Core;
% scaling_factor_protein2_Core;
methods
function MO = poreMaturationReducedPara(N, N1, tpm, ti, force, showplot)
MO@absPoreMaturationMultiPar(N, N1, tpm, ti,force, showplot)
MO.protF = MO.protA;
for i=1:2
MO.protF(i).core = MO.protA(i).tot_core;
MO.protF(i).noncore = MO.protA(i).tot_noncore;
end
MO.parnames = '';
for i=1:2
MO.parnames = [MO.parnames 'kpm' num2str(i) ' '];
end
for i=1:2
MO.parnames = [MO.parnames 'kip' num2str(i) ' '];
end
MO.parnames = [MO.parnames 'kd v fMP_NC fMP_C scal_prot1_NC scal_prot2_NC scal_prot1_C scal_prot2_C'];
MO.LB = [0 0 0 0 0 0 0.9 0.385 0.8 0.8 0.8 0.8];
MO.HB = [10 10 10 10 0.1 0.1 0.94 0.617 1.2 1.2 1.2 1.2 ];
MO.name = mfilename('class');
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)
%% 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
% par(11) - v, production of intermediates
% par(12) - kd, degradation of mature pores
%experimental data has 1 min time resolution
smoothedArea = fnval(MO.protF(1).surfspl, MO.protF(1).time');
parloc = [par(1)*ones(1,MO.N1-1) par(2)*ones(1, MO.N-MO.N1)];
PM = MO.solveModelN([parloc par(5) par(6)], MO.timeI,1);
parlocIP = [par(3)*ones(1,MO.N1-1) par(4)*ones(1, MO.N-MO.N1)];
IP = MO.solveModelN([parlocIP par(5) par(6)], MO.timeI,1);
%compute solution for specific time
PM_t = deval(PM, MO.timeI);
IP_t = deval(IP, MO.timeI);
%% compute protein amounts in non-core and core regions
tend = MO.protF(1).time(end);
iend = MO.protF(1).time(end) + 1;
p1NC_pm = [zeros(1, MO.tpm(1)-4) sum(PM_t(MO.N1:end,1:iend-MO.tpm(1)))]*par(7)*par(9);
p1NC_ip = [zeros(1, MO.ti(1)-4) sum(IP_t(MO.N1:end,1:iend-MO.ti(1)))]*(1-par(7))*par(9);
p1NC = (p1NC_pm+p1NC_ip);
p1NC_den = p1NC./smoothedArea;
p2NC_pm = [zeros(1,MO.tpm(1)-4) PM_t(end,1:iend-MO.tpm(1))]*par(7)*par(10);
p2NC_ip = [zeros(1,MO.ti(1)-4) IP_t(end,1:iend-MO.ti(1))]*(1-par(7))*par(10);
p2NC = (p2NC_pm + p2NC_ip);
p2NC_den = p2NC./smoothedArea;
p1C_pm = [zeros(1,MO.tpm(2)-4) sum(PM_t(MO.N1:end,1:iend-MO.tpm(2)))]*par(8)*par(11);
p1C_ip = [zeros(1,MO.ti(2)-4) , sum(IP_t(MO.N1:end,1:iend-MO.ti(2)))]*(1-par(8))*par(11);
p1C_ip_den = p1C_ip/(1-par(8))/par(11)./smoothedArea*max(smoothedArea);
p1C = (p1C_pm + p1C_ip);
p1C_den = p1C./smoothedArea;
p2C_pm = [zeros(1,MO.tpm(2)-4) PM_t(end,1:iend-MO.tpm(2))]*par(8)*par(12);
p2C_ip = [zeros(1,MO.ti(2)-4) IP_t(end,1:iend-MO.ti(2))]*(1-par(8))*par(12);
p2C_ip_den = p2C_ip/(1-par(8))/par(12)./smoothedArea*max(smoothedArea);
p2C = (p2C_pm+ p2C_ip);
p2C_den = p2C./smoothedArea;
%densities
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), 'C_ip_den', p1C_ip_den);
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), 'C_ip_den', p2C_ip_den);
prot = [p1;p2];
%%
end
function M = modelMatrix(MO, par)