function [ dx ] = bilbico( t, x, k_assEvoTarget, k_dissEvoTarget, k_expressP1, k_expressP2 )
% function for simulating the binding complex
%   input:
%       t: time span
%       x: vector 
%           (1:Evobody (Evo), 2:Target (Target Protein), 3:DNA, 4:Evobody-Target, 5:Target-DNA, 6: Evobody-Target-DNA, 
%            7:bla_Cyt (beta-lactamase in cytoplasm), 8: bla_Peri (beta-lactamase in periplasm), 9:amp_Peri (Ampicillin in periplasm), 
%            10: amp_Medium (Ampicillin in Medium))
%   output:

%definition of parameter
%k1=379000;
%kdiss1= 0.00272;

%association and dissassociation rates for Target and DNA
k_assTargetDNA = 10000;
k_dissTargetDNA = 2*10^-8 * k_assTargetDNA;
%degeneration rates for target and Evobodz
k_degEvo = 4*10^-4;
k_degTarget = 2*10^-3;

%basal expression rate for bla
%k_basalexprbla = 0.01;
k_basalexprbla = 1.84504 * 10^-11;
%expression rate for bla
%k_exprbla = 0.99;
k_exprbla = 1.82659 * 10^-3;
%degeneration rate for bla
k_degbla = 1.82659 * 10^-3;
%transportation rate for bla from cytoplasm into periplasm
k_transportbla = 1*10^-3;

%transport into periplasm
k_transportamp = 1*10^-3;

dx = [0;0;0;0;0;0;0;0;0];

%change of Evobody
%d/dt[Evobody] = k_exprP1[1] + k_dissEvoTarget[Evobody-Target] - k_assEvoTarget[Evobody][Target] + k_dissEvoTarget[Evobody-Target-DNA] - k_assEvoTarget[Evobody][Target-DNA] - k_degEvo[Evobody]
dx(1)=k_expressP1*1+ k_dissEvoTarget*x(4)-k_assEvoTarget*x(1)*x(2)+k_dissEvoTarget * x(6)- k_assEvoTarget*x(1)*x(5)-k_degEvo*x(1);
%change of Target
%d/dt[Target] = k_exprP2[1] + k_dissEvoTarget[Evobody-Target] - k_assEvoTarget[Evobody][Target] + k_dissTargetDNA[Target-DNA] - k_assTargetDNA[Target][DNA] - k_degTarget[Target]
dx(2)= k_expressP2*1+k_dissEvoTarget*x(4)-k_assEvoTarget*x(1)*x(2)+k_dissTargetDNA*x(5)-k_assTargetDNA*x(2)*x(3)-k_degTarget*x(2);
%occupancy of DNA
%d/dt[DNA] = k_dissTargetDNA[Target-DNA] - k_assTargetDNA[Target][DNA] + k_dissTargetDNA[Evobody-Target-DNA] - k_assTargetDNA[Evobody-Target][DNA]
dx(3)= k_dissTargetDNA*x(5)-k_assTargetDNA*x(2)*x(3)+k_dissTargetDNA*x(6)- k_assTargetDNA *x(4)*x(3);
%change of Evobody-Target complex
%d/dt[Evobody-Target] = - k_dissEvoTarget[Evobody-Target] + k_assEvoTarget[Evobody][Target] + k_dissTargetDNA[Evobody-Target-DNA] - k_assTargetDNA[Evobody-Target][DNA]
dx(4)= -k_dissEvoTarget*x(4)+k_assEvoTarget* x(1)*x(2)+k_dissTargetDNA*x(6)-k_assTargetDNA*x(4)*x(3);
%change of Target-DNA complex
%d/dt[Target-DNA] = k_assTargetDNA[Target][DNA] + k_dissEvoTarget[Evobody-Target-DNA] - k_assEvoTarget[Evobody][Target-DNA] - k_dissTargetDNA[Target-DNA]
dx(5)= k_assTargetDNA *x(2)*x(3)+k_dissEvoTarget*x(6)-k_assEvoTarget*x(1)*x(5)-k_dissTargetDNA*x(5);
%change of Evobody-Target-DNA complex
%d/dt[Evobody-Target-DNA] = - k_dissEvoTarget[Evobody-Target-DNA] + k_assEvoTarget[Evobody][Target-DNA] + k_assTargetDNA[Evobody-Target][DNA] - k_dissTargetDNA[Evobody-Target-DNA]
dx(6)= -k_dissEvoTarget*x(6)+k_assEvoTarget*x(1)*x(5)+k_assTargetDNA*x(4)*x(3)-k_dissTargetDNA*x(6) ;
%change of beta-lactamase in cytoplasm
%d/dt[bla_Cyt] = k_exprbla[1] + k_exprblaETD[Evobody-Target-DNA] - k_degbla[bla_Cyt] - k_transpbla[bla_Cyt]
dx(7) = k_basalexprbla*1+k_exprbla*x(6)-k_degbla*x(7)-k_transportbla*x(7);
%change of bla in periplasm
% d/dt[bla_Peri] = - k_degbla[bla_Peri] + k_transpbla[bla_Cyt]

if x(8) < 4.208*10^-6
    dx(8)=-k_degbla*x(8)+k_transportbla*x(7);
else
    dx(8)=0;
end

%change of ampicillin
%   Assosiation rate k_diff (M^-1*s^-1) = 4*pi*r_(A+B)*D'*sin(alpha)*N_A
%   A = molecule A (Ampicillin); B = molecule B (beta-lactamase)
%   r_(A+B) = r_ab = radius A + radius B * 10^(-8)(dm)
%   D' = d_stroke = relative diffusion coefficient (cm^2/s) 
%      = diff coeff A + diff coeff B; must be in dm^2/s -> *10^-2
%   N_A = Avogadro's constant = 6.02214*10^23 mol^-1
%   alpha = angle

%nm -> dm
r_ab = (2.0277+0.4649) * 10^(-8);

%in cm^2/s
d_stroke = (1.237*10^(-8))+(2.825*10^(-9));
%in dm^2/s
d_stroke= d_stroke*10^(-2);

%Avogadro constant
N_A=6.02214*10^(23);

%25 degree
alpha = 5*pi/36;

%change of ampicillin in periplasm
%d/dt[amp_Peri] = - k_diff * [amp_Peri] * [bla] + k_transportamp[amp_Medium] - k_transportamp[amp_Peri]
%dx(9)=-4*pi*((2.0277+0.4649)*10^-8)*(((1.237*10^(-8))+(2.825*10^(-9)))*10^-2)*N_A*x(8)*x(9);
%M^-1*s^-1=l/(mol*s)=dm^3/(mol*s)
kdiff = 4*pi*r_ab*d_stroke*N_A*sin(alpha);
dx(9) = - kdiff * x(8) * x(9)+k_transportamp* 0.000572-k_transportamp*x(9);

%change of ampicillin in medium
% d/dt[amp_Medium] = - k_transportamp[amp_Medium] + k_transportamp[amp_Peri]
dx(10) = -k_transportamp* x(10)+k_transportamp*x(9);

end

