%% VK-model as decribde in van den Herik & ten Tusscher
colors = {'#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd','#8c564b','#e377c2'};

s0 = 0.1:1:2500; %Source/Sink concentration (mM)

%V
v1 = 22.5; v2 = 2.5; k1=75; k2 =75; %V (nmol/s) and K (mM) parameters
u1_v = v1 * s0 ./ (s0 +k1);
u2_v = v2 * s0 ./ (s0 +k2);
u0_v = u1_v + u2_v;
PC1_v = u1_v./u0_v *100;

%K
v1 = 12.5; v2 = 12.5; k1=75; k2 =750; 
u1_k = v1 * s0 ./ (s0 +k1);
u2_k = v2 * s0 ./ (s0 +k2);
u0_k = u1_k + u2_k;
PC1_k = u1_k./u0_k *100;

%R
v1 = 12.5; v2 = 12.5; k1=75; k2 =75; 
u1_r = v1 * s0 ./ (s0 +k1);
u2_r = v2 * s0 ./ (s0 +k2);
u0_r = u1_r + u2_r;
PC1_r = u1_r./u0_r *100;

%% Plot results
figure(1)
%V
subplot(3,3,1)
plot(u0_v,PC1_v,'Linewidth', 1.5, 'color' , colors{4}); hold on
xlabel('Loading rate (nmol/s)'); ylabel('partitioning coefficient (%)')
ylim([45 100]); xlim([0 25])
ax = gca(); ax.FontSize = 9;

subplot(3,3,2)
plot(u0_v,u1_v,'Linewidth', 1.5, 'color' , colors{1}); hold on
plot(u0_v,u2_v,'Linewidth', 1.5, 'color' , colors{2})
xlabel('Loading rate (nmol/s)'); ylabel('Unloading rate (nmol/s)'); xlim([0 25])
ax = gca(); ax.FontSize = 9;

subplot(3,3,3)
plot(s0,u1_v,'Linewidth', 1.5, 'color' , colors{1}); hold on
plot(s0,u2_v,'Linewidth', 1.5, 'color' , colors{2})
xlabel('Sink concentration (mM)'); ylabel('Unloading rate (nmol/s)'); xlim([0 2500])
ax = gca(); ax.FontSize = 9;
%K
subplot(3,3,4)
plot(u0_k, PC1_k, 'Linewidth', 1.5, 'color' , colors{4}); hold on
xlabel('Loading rate (nmol/s)'); ylabel('partitioning coefficient (%)')
ylim([45 100]); xlim([0 25])
ax = gca(); ax.FontSize = 9;

subplot(3,3,5)
plot(u0_k,u1_k,'Linewidth', 1.5, 'color' , colors{1}); hold on
plot(u0_k,u2_k,'Linewidth', 1.5, 'color' , colors{2})
xlabel('Loading rate (nmol/s)'); ylabel('Unloading rate (nmol/s)'); xlim([0 25])
ax = gca(); ax.FontSize = 9;

subplot(3,3,6)
plot(s0,u1_k,'Linewidth', 1.5, 'color' , colors{1}); hold on
plot(s0,u2_k,'Linewidth', 1.5, 'color' , colors{2})
xlabel('Sink concentration (mM)'); ylabel('Unloading rate (nmol/s)'); xlim([0 2500])
ax = gca(); ax.FontSize = 9;
%R
subplot(3,3,7)
plot(u0_r, PC1_r,'Linewidth', 1.5, 'color' , colors{4});; hold on
xlabel('Loading rate (nmol/s)'); ylabel('partitioning coefficient (%)')
ylim([45 100]); xlim([0 25])
ax = gca(); ax.FontSize = 9;

subplot(3,3,8)
plot(u0_r,u1_r, 'Linewidth', 1.5, 'color' , colors{1}); hold on
plot(u0_r,u2_r, 'Linewidth', 1.5, 'color' ,colors{2})
xlabel('Loading rate (nmol/s)'); ylabel('Unloading rate (nmol/s)'); xlim([0 25])
ax = gca(); ax.FontSize = 9;

subplot(3,3,9)
plot(s0,u1_r,'Linewidth', 1.5, 'color' , colors{1}); hold on
plot(s0,u2_r,'Linewidth', 1.5, 'color' , colors{2})
xlabel('Sink concentration (mM)'); ylabel('Unloading rate (nmol/s)'); xlim([0 2500])
ax = gca(); ax.FontSize = 9;
