%% Simulate PET kinetics with 2-tissue compartment model
% Comparison of 6 curves
% Stefanie Krämer, ETH Zurich, 2017

close all;
clearvars;

scatter = 'off'; % chose whether to show some simulated scatter
NumCurves = 6; % Number of curves to show, starting with (1), max 6.
Colors = [0 0 0;0 0 1;0 1 0;1 0 0;1 0 1;0 1 1];

K1(1) = 0.4/60; % black
K1(2) = 2*K1(1); % blue
K1(3) = K1(1); % green
K1(4) = K1(1); % light blue
K1(5) = K1(1); % pink
K1(6) = K1(1); % red

k2(1) = 0.4/60; % black
k2(2) = 2*k2(1); % blue
k2(3) = k2(1); % green
k2(4) = k2(1); % light blue
k2(5) = k2(1); % pink
k2(6) = k2(1); % red

k3(1) = 0.8/60; % black
k3(2) = k3(1); % blue
k3(3) = 2*k3(1); % green
k3(4) = 2*k3(1); % light blue
k3(5) = 0; % pink
k3(6) = k3(1); % red

k4(1) = 0.1/60; % black
k4(2) = k4(1); % blue
k4(3) = k4(1); % green
k4(4) = 2*k4(1); % light blue
k4(5) = k4(1); % pink
k4(6) = k4(1); % red

K1M(1) = 0;
K1M(2) = 0;
K1M(3) = 0;
K1M(4) = 0;
K1M(5) = 0;
K1M(6) = 0.1/60;

k2M(1) = K1M(1);
k2M(2) = 0.01;
k2M(3) = 0.01;
k2M(4) = 0.01;
k2M(5) = 0.01;
k2M(6) = 0.1/60;

Hkt = 0.44;

% Input function Params (kBq/ccm and 1/s)
D(1)=50; % MBq/kg
D(2)=50;
D(3)=50;
D(4)=50;
D(5)=50;
D(6)=50;

tinf(1)=10;
tinf(2)=100;
tinf(3)=200;
tinf(4)=500;
tinf(5)=1000;
tinf(6)=2000;

Va=0.3; % L/kg
Vz=10; % L/kg

CLa(1) = 10000/60/1000/70;
CLa(2) = CLa(1);
CLa(3) = CLa(1);
CLa(4) = CLa(1);
CLa(5) = CLa(1);
CLa(6) = CLa(1);

CLz(1) = 800/60/1000/70; % must be < CLa(1)
CLz(2) = CLz(1);
CLz(3) = CLz(1);
CLz(4) = CLz(1);
CLz(5) = CLz(1);
CLz(6) = CLz(1);

for numSim = 1:numel(K1)
    a(numSim) = CLa(numSim)/Va;
    z(numSim) = CLz(numSim)/Vz;
end

km(1) = 0; % rate constant of parent/total radioactivity in plasma
km(2) = 0;
km(3) = 0;
km(4) = 0;
km(5) = 0;
km(6) = 0.05/60;

fparentinf(1) = 1; % ratio parent/total radioactivity plasma at t = inf.
fparentinf(2) = 1;
fparentinf(3) = 1;
fparentinf(4) = 1;
fparentinf(5) = 1;
fparentinf(6) = 0.2;

v(1)=0.05;
v(2)=v(1);
v(3)=v(1);
v(4)=v(1);
v(5)=v(1);
v(6)=v(1);

% Time (s)
t=(0:1:3600)';
tMin = t/60;
texp=[0 10 20 30 60 90 120 180 240 300 420 540 660 780 900 1200 1500 1800 2100 2400 2700 3000 3300]';

% Plasma function
for numSim = 1:numel(K1)
Plasma(:,numSim) = D(numSim)/tinf(numSim)/Va *...
    ((a(numSim)*Va/Vz-a(numSim))/(a(numSim)*(z(numSim)-a(numSim))) .*...
    (1-exp(-a(numSim).*(t.*(t<=tinf(numSim))+tinf(numSim).*(t>tinf(numSim))))) .*...
    exp(-a(numSim).*(t-tinf(numSim)).*(t>tinf(numSim))) + ...
    (a(numSim)*Va/Vz-z(numSim))/(z(numSim)*(a(numSim)-z(numSim))) .*...
    (1-exp(-z(numSim).*(t.*(t<=tinf(numSim))+tinf(numSim).*(t>tinf(numSim))))) .*...
    exp(-z(numSim).*(t-tinf(numSim)).*(t>tinf(numSim))));
end

% IF of parent and metabolite
for numSim = 1:numel(K1)
    IF(:,numSim) = ...
        Plasma(:,numSim).*((1-fparentinf(numSim)).*exp(-km(numSim).*t) + ...
        fparentinf(numSim));
    IFM(:,numSim) = ...
        Plasma(:,numSim) - IF(:,numSim);
end

% Parent to total radioactivity in plasma
for numSim = 1:numel(K1)
    P2Total(:,numSim) = (1-fparentinf(numSim)).*exp(-km(numSim).*t) + ...
        fparentinf(numSim);
end

% Weighting fuctions
for numSim = 1:numel(K1)
delta1(numSim)=((k2(numSim)+k3(numSim)+k4(numSim))^2-4*k2(numSim)*k4(numSim))^0.5;
theta1(numSim)=(k2(numSim)+k3(numSim)+k4(numSim)+delta1(numSim))/2;
theta2(numSim)=(k2(numSim)+k3(numSim)+k4(numSim)-delta1(numSim))/2;
phi1(numSim)=K1(numSim)*(theta1(numSim)-k3(numSim)-k4(numSim))/delta1(numSim);
phi2(numSim)=K1(numSim)*(theta2(numSim)-k3(numSim)-k4(numSim))/-delta1(numSim);
W(:,numSim)=phi1(numSim)*exp(-theta1(numSim)*t)+phi2(numSim)*exp(-theta2(numSim)*t); % Weighting function
WM(:,numSim) = K1M(numSim)*exp(-k2M(numSim).*t);
end

% Convolution
for numSim = 1:numel(K1)
TAC(:,numSim) =filter(IF(:,numSim),1,W(:,numSim))*(1-v(numSim))+...
    v(numSim).*IF(:,numSim)/(1-Hkt);
TACM(:,numSim) = filter(IFM(:,numSim),1,WM(:,numSim))*(1-v(numSim))+...
    v(numSim).*IFM(:,numSim)/(1-Hkt);
end

% Tissue compartments 1 and 2
for numSim = 1:numel(K1)
TCM2comp1(:,numSim) = (1-v(numSim))*(filter(IF(:,numSim),1,...
    (exp(-(t.*(k2(numSim)+k3(numSim)+k4(numSim)-(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2)...
    .*exp(-(t.*(k2(numSim)+k3(numSim)+k4(numSim)+(k2(numSim)^2+2.*k2(numSim).*k3(numSim)...
    -2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2)...
    .*(K1(numSim).*k3(numSim).*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)-(k2(numSim)^2+2.*k2(numSim).*...
    k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2)...
    +K1(numSim).*k3(numSim).*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)+(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2)...
    +(K1(numSim).*k3(numSim)^2.*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)-(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2))...
    /(k2(numSim)^2+2.*k2(numSim).*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*...
    k4(numSim)+k4(numSim)^2)^(1/2)...
    -(K1(numSim).*k3(numSim)^2.*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)+(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2))...
    /(k2(numSim)^2+2.*k2(numSim).*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+...
    k4(numSim)^2)^(1/2)...
    +(K1(numSim).*k2(numSim).*k3(numSim).*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)-(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2))...
    /(k2(numSim)^2+2.*k2(numSim).*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+...
    k4(numSim)^2)^(1/2)...
    -(K1(numSim).*k3(numSim).*k4(numSim).*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)-(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2))...
    /(k2(numSim)^2+2.*k2(numSim).*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+...
    k4(numSim)^2)^(1/2)...
    -(K1(numSim).*k2(numSim).*k3(numSim).*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)+(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2))...
    /(k2(numSim)^2+2.*k2(numSim).*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+...
    k4(numSim)^2)^(1/2)...
    + (K1(numSim).*k3(numSim).*k4(numSim).*exp((t.*(k2(numSim)+k3(numSim)+k4(numSim)+(k2(numSim)^2+2.*k2(numSim)...
    .*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)+k4(numSim)^2)^(1/2)))/2))...
    /(k2(numSim)^2+2.*k2(numSim).*k3(numSim)-2.*k2(numSim).*k4(numSim)+k3(numSim)^2+2.*k3(numSim).*k4(numSim)...
    +k4(numSim)^2)^(1/2)))/(2.*k3(numSim))));

TCM2comp2(:,numSim) = (1-v(numSim))*TAC(:,numSim) - TCM2comp1(:,numSim);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generate TACexp, TAC1, TAC2 at experimental time points
tIndex=ismember(t,texp);
for numSim = 1:numel(K1)
tTAC=[t TAC(:,numSim)];
texpTAC = tTAC(tIndex,:);
TACexp(:,numSim) = texpTAC(:,2);
end

% Scatter for TACs
Mean=1;
Stdev=0.025;
for numSim = 1:numel(K1)
Scatterfactor(:,numSim) = Stdev.*randn(length(texp),1)+Mean;
ScatterfactorIF(:,numSim) = Stdev.*randn(length(t),1)+Mean;
TACscatter(:,numSim) = Scatterfactor(:,numSim).*TACexp(:,numSim);
IFscatter(:,numSim) = ScatterfactorIF(:,numSim).*IF(:,numSim);
end

% Logan plots
for numSim = 1:numel(K1)
yLogan(:,numSim)=cumtrapz(TAC(:,numSim))./TAC(:,numSim);
xLogan(:,numSim)=cumtrapz(IF(:,numSim))./TAC(:,numSim);
yLoganexp(:,numSim)=cumtrapz(texp,TACscatter(:,numSim))./TACscatter(:,numSim);
xLoganexp(:,numSim)=cumtrapz(texp,IFscatter(tIndex,numSim))./TACscatter(:,numSim);
end

% Logan reference plot, TAC1 is reference region
for numSim = 1:numel(K1)
yLoganRef(:,numSim)=cumtrapz(TAC(:,numSim))./TAC(:,numSim);
xLoganRef(:,numSim)=(cumtrapz(TAC(:,1))+TAC(:,1)/k2(1))./TAC(:,numSim);
yLoganRefexp(:,numSim)=cumtrapz(texp,TACscatter(:,numSim))./TACscatter(:,numSim);
xLoganRefexp(:,numSim)=(cumtrapz(texp,TACscatter(:,1))+TACscatter(:,1)/k2(1))...
    ./TACscatter(:,numSim);
xLoganRefnok2(:,numSim) = cumtrapz(TAC(:,1))./TAC(:,numSim);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Figure W
figure1 = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes1 = axes('Parent',figure1,...
    'Position',[0.15 0.1 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes1,'on');
hold(axes1,'all');
xlabel('Time (min)','FontSize',12,'FontName','Arial');
ylabel('W','FontSize',12,'FontName','Arial');
Ymatrix1 = W(:,1:NumCurves) ; 
plot1 = plot(tMin,Ymatrix1,'Parent',axes1);
for numCurves = 1:NumCurves
set(plot1(numCurves),'LineWidth',2,'LineStyle','-','Color',Colors(numCurves,:));
end
legend1 = legend(plot1);
set(legend1,'Location','northeast','FontSize',12,'EdgeColor',[1 1 1]);

% Figure IF
figure2 = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes2 = axes('Parent',figure2,...
    'Position',[0.15 0.15 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes2,'on');
hold(axes2,'all');
xlabel('Time (min)','FontSize',14,'FontName','Arial');
ylabel('Input function parent (kBq/cm^{3})','FontSize',14,'FontName','Arial');
Ymatrix1 = IF(:,1:NumCurves) ; 
plot1 = plot(tMin,Ymatrix1,'Parent',axes2);
for numCurves = 1:NumCurves
set(plot1(numCurves),'LineWidth',2,'LineStyle','-','Color',Colors(numCurves,:));
end
Ymatrix2=IFscatter(:,1:NumCurves);
plot2 = plot(tMin,Ymatrix2);
for numCurves = 1:NumCurves
set(plot2(numCurves),'MarkerSize',5,'MarkerFaceColor','none','Marker','o',...
    'LineStyle','none','LineWidth',1,'Color',Colors(numCurves,:),...
    'visible',scatter);
end
legend1 = legend(plot1);
set(legend1,'Location','northeast','FontSize',12,'EdgeColor',[1 1 1]);

% Figure Logan
for numCurves = 1:NumCurves
figure4(numCurves) = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes4 = axes('Parent',figure4(numCurves),...
    'Position',[0.15 0.1 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes4,'on');
hold(axes4,'all');
xlabel('(Cumulative plasma (t))/ROI (t))',...
    'FontSize',12,'FontName','Arial');
ylabel('(Cumulative ROI (t))/ROI (t))',...
    'FontSize',12,'FontName','Arial');
plot1 = plot(xLogan(:,numCurves),yLogan(:,numCurves),'Parent',axes4);
set(plot1(1),'LineWidth',2,'LineStyle','-','Color',Colors(numCurves,:));

plot2 = plot(xLoganexp(:,numCurves),yLoganexp(:,numCurves),'Parent',axes4);
set(plot2(1),'MarkerSize',5,'MarkerFaceColor','none','Marker','o',...
    'LineStyle','none','LineWidth',1,'Color',Colors(numCurves,:));
end

% Figure Logan Reference
for numCurves = 1:NumCurves
figure5(numCurves) = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes5 = axes('Parent',figure5(numCurves),...
    'Position',[0.15 0.1 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes5,'on');
hold(axes5,'all');
xlabel('(Cumulative reference (t)+reference/k2) / ROI (t)','FontSize',12,'FontName','Arial');
ylabel('Cumulative ROI (t) / ROI (t))','FontSize',12,'FontName','Arial');
plot1 = plot(xLoganRef(:,numCurves),yLoganRef(:,numCurves),'Parent',axes5);
set(plot1(1),'LineWidth',2,'LineStyle','-','Color',Colors(numCurves,:),...
    'DisplayName','with correction for {\itk}_{2}');
plot2 = plot(xLoganRefexp(:,numCurves),yLoganRefexp(:,numCurves),'Parent',axes5);
set(plot2(1),'MarkerSize',5,'MarkerFaceColor','none','Marker','o',...
    'LineStyle','none','LineWidth',1,'Color',Colors(numCurves,:),...
    'DisplayName','with correction for {\itk}_{2}');
plot3 = plot(xLoganRefnok2(:,numCurves),yLoganRef(:,numCurves),'Parent',axes5);
set(plot3(1),'LineWidth',2,'LineStyle',':','Color',Colors(numCurves,:),...
    'DisplayName','without correction for {\itk}_{2}');
legend1 = legend('show');
set(legend1,'Location','northwest','FontSize',12,'EdgeColor',[1 1 1]);
end


% Figure TAC with radiometabolite
for numCurves = 1:NumCurves
figure7(numCurves) = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes7 = axes('Parent',figure7(numCurves),...
    'Position',[0.15 0.1 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes7,'on');
hold(axes7,'all');
xlabel('Time (min)','FontSize',12,'FontName','Arial');
ylabel('ROI activity (kBq/cm^{3})','FontSize',12,'FontName','Arial');
Ymatrix1=[TAC(:,numCurves) TACM(:,numCurves)]'; 
plot1 = plot(tMin,Ymatrix1,'Parent',axes7);
set(plot1(1),'LineWidth',2,'LineStyle','-','Color',[1 0 0],'DisplayName','Parent');
set(plot1(2),'LineWidth',2,'LineStyle','-','Color',[0 1 0],'DisplayName','Metabolit');
legend1 = legend('show','Location','northeast');
end

% Figure W with radiometabolite
for numCurves = 1:NumCurves
figure8(numCurves) = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes8 = axes('Parent',figure8(numCurves),...
    'Position',[0.15 0.1 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes8,'on');
hold(axes8,'all');
xlabel('Time (min)','FontSize',12,'FontName','Arial');
ylabel('W ','FontSize',12,'FontName','Arial');
Ymatrix1=[W(:,numCurves) WM(:,numCurves)]'; 
plot1 = plot(tMin,Ymatrix1,'Parent',axes8);
set(plot1(1),'LineWidth',2,'LineStyle','-','Color',[1 0 0],...
    'DisplayName','Parent');
set(plot1(2),'LineWidth',2,'LineStyle','-','Color',[0 1 0],...
    'DisplayName','Metabolite');
legend1 = legend(plot1);
set(legend1,'Location','northeast','FontSize',12,'EdgeColor',[1 1 1]);
end

% Figure ratio parent/total in plasma
figure9 = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes9 = axes('Parent',figure9,...
    'Position',[0.15 0.1 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes9,'on');
hold(axes9,'all');
axis([0 60 0 1]);
xlabel('Time (min)','FontSize',12,'FontName','Arial');
ylabel('Parent/total activity in plasma','FontSize',12,'FontName','Arial');
Ymatrix1=P2Total(:,1:NumCurves); 
plot1 = plot(tMin,Ymatrix1,'Parent',axes9);
for numCurves = 1:NumCurves
set(plot1(numCurves),'LineWidth',2,'LineStyle','-','Color',Colors(numCurves,:));
end
legend1 = legend(plot1);
set(legend1,'Location','northeast','FontSize',12,'EdgeColor',[1 1 1]);

% Figure TACs
figure3 = figure('PaperSize',[21 29.7],'Color',[1 1 1]);
axes3 = axes('Parent',figure3,...
    'Position',[0.15 0.15 0.8 0.8],...
    'FontSize',12,...
    'FontName','Arial');
box(axes3,'on');
hold(axes3,'all');
xlabel('Time (min)','FontSize',14,'FontName','Arial');
ylabel('ROI activity (kBq/cm^{3})','FontSize',14,'FontName','Arial');
Ymatrix1 = TAC(:,1:NumCurves); 
plot1 = plot(tMin,Ymatrix1,'Parent',axes3);
for numCurves = 1:NumCurves
set(plot1(numCurves),'LineWidth',2,'LineStyle','-',...
    'Color',Colors(numCurves,:),'DisplayName',num2str(numCurves));
end
Ymatrix2=TACscatter(:,1:NumCurves);
plot2 = plot(texp/60,Ymatrix2,'Parent',axes3);
for numCurves = 1:NumCurves
set(plot2(numCurves),'MarkerSize',5,'MarkerFaceColor','none','Marker','o',...
    'LineStyle','none','LineWidth',1,'Color',Colors(numCurves,:),...
    'visible',scatter);
end
legend1 = legend(plot1);
set(legend1,'Location','northeast','FontSize',12,'EdgeColor',[1 1 1]);

Params = table(K1'*60,k2'*60,k3'*60,k4'*60,K1M'*60,k2M'*60,...
    'VariableNames',{'K1','k2','k3','k4','K1M','k2M'});