pDusk_pDawn_evaluation.m

%% pDusk_pDawn_evaluation
%% Clear History
clc;                % clear display
clear all;          % clear Workspace
close all;          % close all figures 
clear workspace;    % clear workspace


%% 
%% Common Parameters of pDusk and pDawn
time_interval = 0.010;      % define time steps for the figures. 0.01h means taking a time point every 3.6s.
n_pDusk  = 16;              % number of parameters of the ODE system of pDusk
n_pDawn  = 21;              % number of parameters of the ODE system of pDawn
n_parameters = 21;          % n_parameters = n_pDawn (for better readability of code)
n = 10000;                   % number of the sample parameter sets to be generated
System = 'pDusk_pDawn';     % for writing data files
Underscore = '_';           % for writing data files


%% Reading the Generated Data
Underscore = '_';

% Output_pDusk
spec = '%d%s%.3f%s%s%sOutput_pDusk.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
Output_pDusk = csvread(str);

% Output_pDawn
spec = '%d%s%.3f%s%s%sOutput_pDawn.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
Output_pDawn = csvread(str);

% BestOutput_pDusk
spec = '%d%s%.3f%s%s%sBestOutput_pDusk.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
BestOutput_pDusk = csvread(str);

% BestOutput_pDawn
spec = '%d%s%.3f%s%s%sBestOutput_pDawn.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
BestOutput_pDawn = csvread(str);

% Score
spec = '%d%s%.3f%s%s%sScore.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
Score = csvread(str);

% Parameter_Set
spec = '%d%s%.3f%s%s%sParameter_Set.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
Parameter_Set = csvread(str);

n = 1000;
% pDusk Ratio
spec = '%d%s%.3f%s%s%spDusk_Ratio.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
pDusk_Ratio = csvread(str);

% pDawn Ratio
spec = '%d%s%.3f%s%s%spDawn_Ratio.csv'; 
str = sprintf(spec, n, Underscore, time_interval, Underscore, System, Underscore);
pDawn_Ratio = csvread(str);
n = 10000;


%% Unit Conversion of Light Intensities
% As the ODE system follows the units µmol, meters (m) and hours (h), the
% given light intensities from Ohlendorf et al. (2012) in [µW cm^-2] needed
% to be transformed. This was done according to: Berthold Technologies 
% (https://www.berthold.com/en/bio/how-do-i-convert-irradiance-photon-flux)
% Addition: There is a correction to the given Ohlendorf nW cm^-2 given. 
% It should read µW cm^-2 (Corrigendum: DOI: 10.1016/j.jmb.2013.11.005).

% Constants
h = 6.63*10^-34;                        % [Js] Planck's constant
c = 2.988*10^8;                         % [m/s] speed of light
lambda_nm = 470;                        % [nm] wavelength
N_A = 6.022*10^23;                      % [mol^-1] Avogadro's constant


%% 
%% pDusk
% Defining variables for pDusk
t   = 0:time_interval:17;   % Ohlendorf et al. (2013) induced the bacteria for 17h, i.e. 61200s. Thus the same time frame is used
x0  = [1 0 0 0 1 0 1];      % initial guesses for the variables of the ODE system. In the beginning state, the system is assumed to be in dark conditions, thus y_DD = 100% and j_a = 100% ("all" inactive FixJ (j_i) is phosphorylated)
p = zeros(n_parameters);    % creating a zero-filled parameter vector


%% Loading the Data Points of Ohlendorf et al. (2012) for pDusk
% The data points of figure 2a in Ohlendorf et al. (2012) were retrieved
% through WebPlotDigitizer (http://arohatgi.info/WebPlotDigitizer) by (Ankit
% Rohatgi) to be used for the scoring function to evaluate the parameter
% sets generated by the latin hypercube sampling.
DataPoints_pDusk = csvread('DataPoints_pDusk.csv');     % reading the data points of figure 2a from Ohlendorf et al. (2012)
lightintensities_pDusk = [DataPoints_pDusk(:,1)'];      
RFP_datapoints_pDusk = [DataPoints_pDusk(:,2)'];

% Conversion of Light Intensities
N_pDusk = lightintensities_pDusk*10^-2;                 % µW*cm^-2 --- *10^-2 ---> W*m^-2
E_p = h*c/lambda_nm;                                    % distinct energy quanta of a photon
N_p = N_pDusk*lambda_nm*1E6/(1.988*1E-16);              % number of photons; lambda in [nm]
E_q = N_p/N_A;                                      % photonflux N = E_q in [µmol*m^-2*s^-1]
E_q_h = E_q*60*60;                                      % photonflux in [µmol*m^-2*h^-1]
N_pDusk = E_q_h;


%% Loading the Standard Deviation Values from Ohlendorf et al. (2012)
F = csvread('Standard_Deviation_Ohlendorf_pDusk.csv');
StDev_Ohlendorf_pDusk = [F(:,2)'];                                                            % this is StDev in logscale
StDev_Sorted_pDusk = StDev_Ohlendorf_pDusk([1,3,5,7,9,11,13,15]);                             % subsetting upper bound StDev values
StDev_Sorted_pDusk = [StDev_Sorted_pDusk; StDev_Ohlendorf_pDusk([2,4,6,8,10,12,14,16])];      % subsetting lower bound StDev values
% StDev_RFP_pDusk = (StDev_Sorted_pDusk(1,:)-StDev_Sorted_pDusk(2,:))/2;                      % extracting StDev values (not considering log scale)
StDev_RFP_pDusk = (exp(StDev_Sorted_pDusk(1,:))-exp(StDev_Sorted_pDusk(2,:)))/2;              % extracting StDev values (considering log scale)


%%
%% pDawn
% Defining variables for pDawn
t   = 0:time_interval:17;   % Ohlendorf et al. (2013) induced the bacteria for 17h, i.e. 61200s. Thus the same time frame is used
x0  = [1 0 0 0 1 0 0 0 0];  % initial guesses for the variables of the ODE system. In the beginning state, the system is assumed to be in dark conditions, thus y_DD = 100% and j_a = 100% ("all" inactive FixJ (j_i) is phosphorylated)
p = zeros(n_parameters);    % creating a zero-filled parameter vector


%% Loading the Data Points of Ohlendorf et al. (2012) for pDawn
% The data points of figure 2a in Ohlendorf et al. (2012) were retrieved
% through WebPlotDigitizer (http://arohatgi.info/WebPlotDigitizer) by (Ankit
% Rohatgi) to be used for the scoring function to evaluate the parameter
% sets generated by the latin hypercube sampling.
DataPoints_pDawn = csvread('DataPoints_pDawn.csv');     % reading the data points of figure 2a from Ohlendorf et al. (2012)
lightintensities_pDawn = [DataPoints_pDawn(:,1)'];      
RFP_datapoints_pDawn = [DataPoints_pDawn(:,2)'];

% Conversion of Light Intensities
N_pDawn = lightintensities_pDawn*10^-2;                 % µW*cm^-2 --- *10^-2 ---> W*m^-2
E_p = h*c/lambda_nm;                                    % distinct energy quanta of a photon
N_p = N_pDawn*lambda_nm*1E6/(1.988*1E-16);              % number of photons; lambda in [nm]
E_q = N_p/N_A;                                      % photonflux N = E_q in [µmol*m^-2*s^-1]
E_q_h = E_q*60*60;                                      % photonflux in [µmol*m^-2*h^-1]
N_pDawn = E_q_h;


%% Loading the Standard Deviation Values from Ohlendorf et al. (2012)
F = csvread('Standard_Deviation_Ohlendorf_pDawn.csv');
StDev_Ohlendorf_pDawn = [F(:,2)'];                                                            % this is StDev in logscale
StDev_Sorted_pDawn = StDev_Ohlendorf_pDawn([1,3,5,7,9,11,13,15]);                             % subsetting upper bound StDev values
StDev_Sorted_pDawn = [StDev_Sorted_pDawn; StDev_Ohlendorf_pDawn([2,4,6,8,10,12,14,16])];      % subsetting lower bound StDev values
StDev_RFP_pDawn = (exp(StDev_Sorted_pDawn(1,:))-exp(StDev_Sorted_pDawn(2,:)))/2;              % extracting StDev values (considering log scale)


%%
%% Common Analysis
%% Selecting the Parameters of the Best Parameter Set
BestScore = Score(Score(:,2) == min(Score(:,2)),:);
[val,ord]=sort(Score(:,2));
sortedScores=Score(ord,:);
best_p = Parameter_Set(sortedScores(1,1),:);
best_p_pDusk = best_p(1:16);
best_p_pDawn = best_p(1:end);


%% creating subsets
% pDusk
constrained_pDusk_Ratio_bigger_one = pDusk_Ratio(pDusk_Ratio(:,3) > 1 & pDusk_Ratio(:,2) == 1,:);    % A/T > 1       at N(1)
constrained_pDusk_Ratio_smaller_one = pDusk_Ratio(pDusk_Ratio(:,3) < 1 & pDusk_Ratio(:,2) == 8,:);   % A/T < 1       at N(8)
subset = ismember(constrained_pDusk_Ratio_bigger_one(:,1), constrained_pDusk_Ratio_smaller_one(:,1));
subset = constrained_pDusk_Ratio_bigger_one(subset,1)';


%% 
%% Separate Analysis
%% Simulating Each Parameterset (j) per Each Light Intensity N(i)
% Defining variables for pDusk
t   = 0:time_interval:17;   % Ohlendorf et al. (2013) induced the bacteria for 17h, i.e. 61200s. Thus the same time frame is used
x0  = [1 0 0 0 1 0 0];      % initial guesses for the variables of the ODE system. In the beginning state, the system is assumed to be in dark conditions, thus y_DD = 100% and j_a = 100% ("all" inactive FixJ (j_i) is phosphorylated)
p = zeros(n_parameters);    % creating a zero-filled parameter vector

% best parameter set
Xvalues = [];
Ni = [];
Time = [];
N_pDusk_0 = [0 N_pDusk];
for i = 1:length(N_pDusk_0)
    [T, X] = ode45(@(t,x) pDusk_function(t,x,best_p_pDusk,N_pDusk_0(i)), t, x0);
    Xvalues = [Xvalues; X];
    
    J = 1:length(X(:,1));
    J(J~=0) = i;
    Ni = [Ni; J'];
    
	Time = [Time; T];
end 
X = [];
X_best_p_pDusk = [Ni, Time, Xvalues(:,1:7)]; 


%% Plotting the response of all system variables of pDusk
figure(1);
for i = 1:length(N_pDusk);
    % plotting the mean response with 95% "confidence interval"
    subplot(3,3,i);
    % best parameter set
    red = plot(T,X_best_p_pDusk(X_best_p_pDusk(:,1) == i, end) );
    str = sprintf('N = %f',N_pDusk(i));
    str = {'Response of Species in pDusk per Light Intensity', str};
    title(str);
    xlabel('time [h^1]');
    ylabel('concentration [µmol]');
end

% adding a legend to the graph
subplot(3,3,9);
for i = 1;
    white = plot(T(1),X_best_p_pDusk(1, end));
    % adding Legend
    str = {'Legend'};
    title(str);
    legend ('RFP_p'); 
end

print('-painters', '-dpng', '-f1');
%% pDawn
% Defining variables for pDawn
t   = 0:time_interval:17;   % Ohlendorf et al. (2013) induced the bacteria for 17h, i.e. 61200s. Thus the same time frame is used
x0  = [1 0 0 0 1 0 0 0 0];  % initial guesses for the variables of the ODE system. In the beginning state, the system is assumed to be in dark conditions, thus y_DD = 100% and j_a = 100% ("all" inactive FixJ (j_i) is phosphorylated)
p = zeros(n_parameters);    % creating a zero-filled parameter vector

% Simulating Each Parameterset (j) per Each Light Intensity N(i) pDawn
% output = Output_pDawn(Output_pDawn(:,1) == sortedScores(1,1), end-2:end-1);
% best parameter set
Xvalues = [];
Ni = [];
Time = [];
N_pDawn_0 = [0 N_pDawn];
for i = 1:length(N_pDawn_0)
%     [T, X] = ode45(@(t,x) pDawn_function(t,x,best_p_pDawn,N_pDawn(i)), t, x0);
    [T, X] = ode45(@(t,x) pDawn_function(t,x,best_p_pDawn,N_pDawn_0(i)), t, x0);
    Xvalues = [Xvalues; X];
    
    J = 1:length(X(:,1));
    J(J~=0) = i;
    Ni = [Ni; J'];
    
	Time = [Time; T];
end 
X = [];
X_best_p_pDawn = [Ni, Time, Xvalues(:,1:9)]; 


%% pDusk response for N = 0
Xvalues = [];
Ni = [];
Time = [];
    [T, X] = ode45(@(t,x) pDusk_function(t,x,best_p_pDusk,420), t, x0);
    Xvalues = [Xvalues; X];
    
    J = 1:length(X(:,1));
    J(J~=0) = i;
    Ni = [Ni; J'];
    
	Time = [Time; T];
X = [];
X_best_p_pDusk_N0 = [Ni, Time, Xvalues(:,1:7)];


%% pDawn response for N = 0
Xvalues = [];
Ni = [];
Time = [];
%     [T, X] = ode45(@(t,x) pDawn_function(t,x,best_p_pDawn,N_pDawn(i)), t, x0);
    [T, X] = ode45(@(t,x) pDawn_function(t,x,best_p_pDawn,420), t, x0);
    Xvalues = [Xvalues; X];
    
    J = 1:length(X(:,1));
    J(J~=0) = i;
    Ni = [Ni; J'];
    
	Time = [Time; T];
X = [];
X_best_p_pDawn_N0 = [Ni, Time, Xvalues(:,1:9)];


%% Plotting the response of all system variables of pDawn
figure(2);
for i = 1:length(N_pDawn);
    % plotting the mean response with 95% "confidence interval"
    subplot(3,3,i);
    % best parameter set
    red = plot(T,X_best_p_pDawn(X_best_p_pDawn(:,1) == i, end) );
    clear ylim;
    ylim([0 0.06]);
    str = sprintf('N = %f',N_pDawn(i));
    str = {'RFP Response in pDawn', 'per Light Intensity', str};
    title(str, 'FontSize', 14);
    xlabel('Time [h]', 'FontSize', 12);
    ylabel('Concentration [µmol]', 'FontSize', 12);
    set(findall(gca, 'Type', 'Line'),'LineWidth',1);
    res = 300;
end

print('-painters', '-dpng', '-f2');
%% T--Wageningen_UR--pDusk_pDawn_insilico_experiment
figure(3)
i = 3;
subplot(2,2,1);
plot(T,X_best_p_pDusk_N0(:, end) );
% plot(T,X_best_p_pDusk(X_best_p_pDusk(:,1) == i, end) );
clear ylim;
ylim([0 0.06]);
str = {'RFP Response in pDusk at', 'N = 0 [µmol*m^-^2*h^-^1]'};
title(str, 'FontSize', 14);
xlabel('Time [h]', 'FontSize', 12);
ylabel ('Concentration [µmol]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;

i = 4;
subplot(2,2,2);
plot(T,X_best_p_pDawn_N0(:, end));
% plot(T,X_best_p_pDawn(X_best_p_pDawn(:,1) == i, end) );
clear ylim;
ylim([0 0.06]);
str = {'RFP Response in pDawn at', 'N = 0 [µmol*m^-^2*h^-^1]'};
title(str, 'FontSize', 14);
xlabel('Time [h]', 'FontSize', 12);
ylabel('Concentration [µmol]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;

i = 8;
subplot(2,2,3);
plot(T,X_best_p_pDusk(X_best_p_pDusk(:,1) == i, end) );
clear ylim;
ylim([0 0.06]);
str = sprintf('N = %f [µmol*m^-^2*h^-^1]',N_pDusk(i));
str = {'RFP Response in pDusk at', str};
title(str, 'FontSize', 14);
xlabel('Time [h]', 'FontSize', 12);
ylabel ('Concentration [µmol]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;

i = 8;
subplot(2,2,4);
plot(T,X_best_p_pDawn(X_best_p_pDawn(:,1) == i, end) );
clear ylim;
ylim([0 0.06]);
str = sprintf('N = %f [µmol*m^-^2*h^-^1]',N_pDawn(i));
str = {'RFP Response in pDawn at', str};
title(str, 'FontSize', 14);
xlabel('Time [h]', 'FontSize', 12);
ylabel('Concentration [µmol]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;

print('-painters', '-dpdf', '-f3');
%%
%% Commmon Analysis/ Parameters
%% Comparison 17h lab experiment
pDusk_17 = X_best_p_pDusk(X_best_p_pDusk(:,2) == 17, end);
pDawn_17 = X_best_p_pDawn(X_best_p_pDawn(:,2) == 17, end);

figure(4)
plot(N_pDusk_0, pDusk_17)
hold on
plot(N_pDawn_0, pDawn_17)
hold off
clear ylim
ylim([0 0.06])
str = ('per Light Intensity N');
str = {'RFP Response in pDusk and pDawn'; str};
title(str, 'FontSize', 14);
xlabel('Light Intensities N [µmol*m^-2*h^-1]', 'FontSize', 12);
ylabel('Concentration [µmol]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;
legend('pDusk', 'pDawn');

print('-painters', '-dpng', '-f4');
%% T--Wageningen_UR--pDusk_pDawn_Match_to_Ohlendorf
% Control Plot pDusk_pDawn
figure(5)
% pDusk
RFP_p_computed_pDusk = Output_pDusk(Output_pDusk(:,1) == BestScore(1), end-1); % taking the normalized concentrations of RFP_p depending on the ligth intensities of the parameter set with the best Score_pDusk
subplot(1,2,1);
% Plotting the Ohlendorf DataPoints for pDusk
errorbar(N_pDusk, RFP_datapoints_pDusk, StDev_RFP_pDusk,'.k','LineWidth',1, 'MarkerSize', 10);
set(gca,'YScale','log')
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;
hold on
% Plotting the Computed RFP_p Levels for pDusk
scatter(N_pDusk, RFP_p_computed_pDusk)
plot(N_pDusk, RFP_p_computed_pDusk)
clear ylim;
ylim([0.0008 1.15]);
clear xlim;
xlim([-300 15000]);
title('pDusk', 'FontSize', 14);
xlabel('Light Intensity N in [µmol*m^-^2*h^-^1]', 'FontSize', 12);
ylabel('FL/OD600 [a.u.]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;
hold off

% pDawn
RFP_p_computed_pDawn = Output_pDawn(Output_pDawn(:,1) == BestScore(1), end-1); % taking the normalized concentrations of RFP_p depending on the ligth intensities of the parameter set with the best Score_pDawn
subplot(1,2,2);
% Plotting the Ohlendorf DataPoints for pDawn
errorbar(N_pDawn, RFP_datapoints_pDawn, StDev_RFP_pDawn,'.k','LineWidth',1, 'MarkerSize', 10);
set(gca,'YScale','log')
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;
hold on
% Plotting the Computed RFP_p Levels for pDawn
scatter(N_pDawn, RFP_p_computed_pDawn)
plot(N_pDawn, RFP_p_computed_pDawn)
clear ylim;
ylim([0.0008 1.15]);
clear xlim;
xlim([-300 15000]);
title('pDawn', 'FontSize', 14);
xlabel('Light Intensity N in [µmol*m^-^2*h^-^1]', 'FontSize', 12);
ylabel('FL/OD600 [a.u.]', 'FontSize', 12);
set(findall(gca, 'Type', 'Line'),'LineWidth',1);
res = 300;
hold off

print('-painters', '-dpng', '-f5');