function VirusInfectionFinal
clc
close all
clear all

% definition of a characteristic time scale
tau = 1000;
% definition of a characteristic rate constant
p = 1/tau;
% initial extracellular virus concentration
V0 = 10E7;
% initial receptor concentration
R0 = 10E6;

tic
figure('Name','VirusInfection')
VirusInfectionLMA(tau,V0,R0,p)
figure('Name','VirusInfectionLMA')
VirusInfectionModLMA(tau,V0,R0,p)
toc

function VirusInfectionLMA(tau,V0,R0,p)
% Definition of the rate constants vector
k=zeros(7,2);
% set all rate constants of the forward reactions
k(:,1)=30*p;
% set all rate constants of the reverse reactions
k(:,2)=3*p;
% defining the rate constants for receptor virus complex formation
k(1,1) = 0.5*k(1,1);
% defining the rate constants for receptor unbinding
k(1,2) = 0.5*k(1,2);
% defining the rate constants for receptor dimer virus complex formation
k(2,1) = 1.5*k(2,1);
% defining the rate constants for receptor dimer unbinding
k(2,2) = 1.5*k(2,2);
% rate constants for protein synthetization
k(5,1) = k(5,1)*0.5;
% DNA degradation
k(6,1) = k(6,1)*0.3;
% Virus elimination
k(7,1) = k(7,1)*0.7;

% initial concentrations
y0 = zeros(7,1);
y0(1,1) = V0;         %V_ext0
y0(2,1) = R0;          %R0

% modification degree
m = -1; % (modification not taken into account in this part of the model)
% integration time
tspan = [0 tau];

% integration of odes
[t,y] = ode15s(@odeLMA,tspan,y0,[],k,m,p);

% total receptor concentration
R_sum = y(:,2)+y(:,3)+2*y(:,4);

plot(t,y,'o-')
ylabel('concentration [1/ ml]');
xlabel('time [min]');
ylim([-1 2*10E6])
legend('V_{ext}','R','VR','VR_{2}','V_{nuc}','ssDNA_{nuc}','E_{cyt}')


function VirusInfectionModLMA(tau,V0,R0,p)

% Definition of the rate constants vector
k=zeros(7,2);
% set all rate constants of the forward reactions
k(:,1)=30*p;
% set all rate constants of the reverse reactions
k(:,2)=3*p;
% defining the rate constants for receptor virus complex formation
k(1,1) = 0.5*k(1,1);
% defining the rate constants for receptor unbinding
k(1,2) = 0.5*k(1,2);
% defining the rate constants for receptor dimer virus complex formation
k(2,1) = 1.5*k(2,1);
% defining the rate constants for receptor dimer unbinding
k(2,2) = 1.5*k(2,2);
% rate constants for protein synthetization
k(5,1) = k(5,1)*0.5;
% DNA degradation
k(6,1) = k(6,1)*0.3;
% Virus elimination
k(7,1) = k(7,1)*0.7;
k
% initial concentrations
y0 = zeros(7,1);
y0(1,1) = V0;         %V_ext0
y0(2,1) = R0;         %R0

% modification degree
m = 0:0.05:1;
% integration time
tspan = [0 2*tau];

% initialization of yend vector
yend = zeros(length(m),1);

% loop to calculate the steady state enzyme concentration (E_cyt) for different
% degrees of modification (m)
for i=1:length(m)
    % integration of ODEs
    [t,y] = ode15s(@odeLMA,tspan,y0,[],k,m(i),p);
    subplot(2,1,1)
    plot(t,y(:,7),'--k')
    ylabel('Enzyme concentration');
    xlabel('time [min]');
    legend('E_{cyt}')
    hold on
    
    yend(i,1) = y(end,7);
end
hold off

% plot the steady state enzyme concentration as a function of m
subplot(2,1,2)
plot(m*100,yend,'--k')
ylabel('Enzyme concentration');
xlabel('modification degree [%]');
hold on

% calculate the optimal modification degree
[~,j] = max(yend);
optimum = [m(j)*100 yend(j)]
% plot the optimum
plot(m(j)*100,yend(j),'or','MarkerSize',10)
hold off

function dydt = odeLMA(~,y,k,m,p)
% extracellular virus concentration
V_ext       =   y(1);
% Receptor concentration
R           =   y(2);
% receptor monomer virus complex
VR          =   y(3);
% receptor dimer virus complex
VR_2        =   y(4);
% intranuclear virus concentration
V_nuc       =   y(5);
% intranuclear ssDNA concentration
ssDNA_nuc   =   y(6);
% intracellular enzyme concentration
E_cyt       =   y(7);

% initialization of dydt
dydt = zeros(7,1);
% dV_ext/dt
dydt(1,1) = 	-k(1,1)*V_ext*R +k(1,2)*VR -k(7,1)*V_ext;
% dR/dt
dydt(2,1) = 	-k(1,1)*V_ext*R +k(1,2)*VR -k(2,1)*VR*R +k(2,2)*VR_2 +k(3,1)*kmod(m)*2*VR_2;
% dVR/dt
dydt(3,1) = 	-k(1,2)*VR +k(1,1)*V_ext*R -k(2,1)*VR*R +k(2,2)*VR_2;
% dVR_2/dt
dydt(4,1) =  	-k(2,2)*VR_2 +k(2,1)*VR*R -k(3,1)*VR_2;
% dV_nuc/dt
dydt(5,1) = 	-k(4,1)*V_nuc +k(4,2)*ssDNA_nuc +k(3,1)*kmod(m)*VR_2;
% dssDNA_nuc/dt
dydt(6,1) = 	-k(4,2)*ssDNA_nuc +k(4,1)*V_nuc*rmod(m) -k(5,1)*ssDNA_nuc -k(6,1)*ssDNA_nuc;
% dE_cyt/dt
dydt(7,1) = 	k(5,1)*ssDNA_nuc;

function k=kmod(m)
% rate consant as a function of m
l=0.5;
a=0.1;
if m>=a
    k = (m-a)^2/(l^2+(m-a)^2);
elseif m==-1
    k = 1;
else
    k = 0;
end

function r=rmod(m)
% fraction of filled capsids as a function of m
l=0.3;
a=0.15;
if m>=a
    r = 1-(m-a)^2/(l^2+(m-a)^2);
else
    r = 1;
end

