% Conor O'Brien, Qi Zhang, Prodromos Daoutidis, and Wei-Shou Hu

% This script is used to simulate the in-silico process model

clear
close all
clc

% ------------------------------------------------------------------------
%             Experimental data entry/import
% ------------------------------------------------------------------------

% Example starting conditions for simulation
% [Glucose, lactate, VCD, osmolarity]
% Row 1: low titer run
% Row 2: high titer run
exptStartingConditions = [
    58.4, 17.0, 1.2, 316.4
    58.2, 14.4, 1.2, 315.9];
index = [1 100]; % These numbers should correspond to the rank of each run in the list of starting conditions. Length must be the same as the number of rows of experimental starting points
totalLength = 100; % This number corresponds to the total number of runs in the dataset
time = [0 260]; % Time in hours over which culture is simulated


% ------------------------------------------------------------------------
%             Initialize model inputs
% ------------------------------------------------------------------------

% Set the time for the initial feed
Tfeed = 90; % Time (h)
% Amount to feed for the timed glucose feed
Tglucose = 22;

% Glucose concentration below which feeding occurs
Gfeed = 15; % [mM]
% Amount to feed for the low glucose feed
Gglucose = 10;


% Fitted parameters for growth and metabolism
% [mumax , kdmax , Kelac , CIexp , qI, Kos, Cos, cellosmo, kla]
growthKinetics = [0.0312814250209082,0.00981732643944096,6057.95203717527,1.42251654503164,0.00323563622939188,0.00229643099389660,286.420442852993,10.56,132];
% [qcell, Kampk, Kstress, E0perm , E0hk, E0pgi, E0pfkfb, E0pfk, E0ald, E0tpi, E0gapd, E0pgk, E0pgm, E0en, E0pk, E0ldh, E0mct, E0g6pd, E0pdhc, Ncd, KBP]
metabolismKinetics = [0.00193102167846227,334.005372644765,2.15895897600695,1,0.306258642306509,3.03209128349985,1,26.0246271883227,2.88731895534180,0.254882322954437,2.43265778097666,1,1,1,0.199008259733499,1.96210038707157,1,1,1,0.48];
% [low KBP, high KBP]
KBPrange = [98.88,213.84];

% Set the initial KBP based on the relative rank of the runs corresponding
% to the starting conditions
KBPculture = rescale(1-((index-1)./(totalLength-1)), KBPrange(1), KBPrange(2),'InputMin',0,'InputMax',1);


% ------------------------------------------------------------------------
%             Culture simulation
% ------------------------------------------------------------------------

% Simulate cultures in order using assigned KBP
C = {};
t = {};

% Function which sets the model to time out after a certain period of time
% has elpased, and to stop evaluation for feeding as needed
timeOverFcn = @(T, Y) OdeTimer(T, Y, Tfeed, Gfeed);
odeopts = odeset('Events',timeOverFcn);

% Simulate cultures
for i = 1:length(index)
    metabolismKinetics(21) = KBPculture(i);
    for j = 1:length(exptStartingConditions(:,1))
        exptToUse = exptStartingConditions(i,:);
    end
    C0 = primeC(growthKinetics,metabolismKinetics,odeopts);
    [ti,Ci] = cultureSim(growthKinetics,metabolismKinetics,exptToUse,time,C0,odeopts, Tglucose, Gglucose);
    t = [t,ti];
    C = [C,Ci];
end


% ------------------------------------------------------------------------
%             Visualization of model simulation
% ------------------------------------------------------------------------

close all

% Color for figure lines based on index
clabelAll = [linspace(1,0,totalLength)', linspace(0,0,totalLength)', linspace(0,1,totalLength)'];
clabel = clabelAll(index,:);

% Generate plots
plotlist = {'Glucose','Lactate','VCD','Osmolarity'};
axislabel = {'Glucose [mM]','Lactate [mM]','VCD [10^6 cells/mL]','Osmolarity [mOsm]'};
modelIndex = [1 15 43 45];
% Set range for the plots, overwrite if changing model parameters
% substantially
plotRange = [0 70 % Glucose
             0 140 % Lactate
             0 8 % VCD
             280 480]; % Osmolarity
fig = figure
for i = 1:length(plotlist)
    a = subplot(2,length(plotlist),i+4)
    set(a,'TickDir','out','Position',[0.23*(i-1)+0.05,0.15,0.18,0.75]);
    modelPlot(t,C,modelIndex(i),clabel,plotlist{i},axislabel{i},plotRange(i,:));
end
set(fig,'Units','normalized','Position',[0.3 0.3 0.45 0.2],'renderer','Painters')


% ------------------------------------------------------------------------
%             Functions
% ------------------------------------------------------------------------

% Prime the starting reactor concentrations to prevent fast transients
function C = primeC(growthKinetics,metabolismKinetics,options)
    Ceglc0 = 58;
    Celac0 = 15;
    Ceosm = 315.7;
    C0=[Ceglc0,0.13638,0.46718,0.063209,10.207,0.070862,0.018154,0.00037755,0.1969,0.033378,0.093299,0.15566,0.07797,33.483,Celac0,0.029469,0.010086,0.00072961,0.003727,0.0068285,0.024475,0.010286,0.0053676,0.00028783,0.15698,6.369,0.30137,0.14727,0.098694,0.2343,1.0581,0.012215,0.058097,47.564,0.00024719,0.40696,0.098211,2.4297,9.7931,5.7984,2.1203,0.012283,1,0,Ceosm,0.12,0.003187,4.78E-06,1.26E-07,0.01,3.92E-06,0];
    tic
    [~,Ctmp] = ode15s(@(t,C) bioprocessmodel(t,C,0,growthKinetics,metabolismKinetics),[0 20],C0,options);
    C = Ctmp(end,:);
end

% Simulate reactor for a set of starting points
function [t,c] = cultureSim(growthKinetics,metabolismKinetics,exptDat,time,C0i,options, Tglucose, Gglucose)
    % Take initial conditions and overwrite with matching experimental
    % values
    C0 = C0i;
    C0(1) = exptDat(1);
    C0(15) = exptDat(2);
    C0(43) = exptDat(3);
    C0(45) = exptDat(4);

    tic;

    curT = time(1);
    t = [];
    c = [];
    % Each bioprocess simulation will take multiple integrations, as it the
    % integration is stopped and restarted with every feed during the
    % culture
    for j = 1:15
        [T,C,~,~,ie] = ode15s(@(t,C) bioprocessmodel(t,C,0,growthKinetics,metabolismKinetics),[curT time(end)],C0,options); % Different met kinetics for each iteraction

        % Add end concentrations to vector, prepare to reinitialize
        C0 = C(end,:);
        t = [t;T];
        c = [c;C];

        if ie == 1 % Integration failure, abort here
            break
        elseif ie == 2 % Need glucose addition, restart with additional glucose
            curT = T(end) + 0.001;
            C0(1) = C0(1) + Gglucose;
            C0(45) = C0(45) + Gglucose; % Adjusting osmo
        elseif ie == 3 % Scheduled feed, add glucose
            curT = T(end) + 0.001;
            C0(1) = C0(1) + Tglucose;
            C0(45) = C0(45) + Tglucose; % Adjusting osmo
        elseif T(end) >= time(end)
            break
        end
    end
end

% This function is used both to stop the simulation if it is time for a scheduled feed, 
% if glucose concentration drops too low (triggering another feed), or becomes difficult
% to solve
function [value, isterminal, direction] = OdeTimer(T, Y, Tfeed, Gfeed)
    TimeOut = 10;
    % Tfeed = 70; % change this to function input
    % [(1) cut off evaluation if too long
    %  (2) if glucose too low, stop evaluate to feed glucose
    %  (3) feed glucose at specified time]
    value = [toc-TimeOut, Y(1)-Gfeed, Tfeed-T];
    isterminal = [1 1 1];
    direction = [0 0 0];
end

% Function to produce row of plots for modeling data
function [fig] = modelPlot(t,C,modelIndex,clabel,plotLabel,axislabel,plotRange)
    % Draw one line for each run
    for j = 1:length(t)
        fig = plot(t{j},C{j}(:,modelIndex),'Color',clabel(j,:));
        hold on
    end
    
    title(['Model Fit ' plotLabel])
    
    a = gca;
    % set box property to off and remove background color
    set(a,'box','off','color','none')
    % create new, empty axes with box but without ticks
    b = axes('Position',get(a,'Position'),'box','on','xtick',[],'ytick',[]);
    % set original axes as active
    axes(a)
    % link axes in case of zooming
    linkaxes([a b])
    a.TickDir = 'out'
    
    xlim([0 260])
    if length(plotRange) == 2
        ylim(plotRange)
    else
        ylim auto
    end
    xlabel('Time [hr]')
    ylabel(axislabel)
end