function [dx_dt]=final_setup_dsf(t,x,mu,clp)

% Gillian Forsyth, Aleksandra Plochocka, Robyn Shuttleworth 
%
% The following function, final_setup_dsf, calls the function,
% final_solver_dsf, in order to solve the deterministic system of equations
% modelling the DSF pathway our engineered E.coli cell. The pathway is described by the
% following equations:

% RpfG > RpfG[P] at a rate k1*D
% RpfG[P] + c-di-GMP > RpfG at a rate k2
% Clp + c-di-GMP > Clp.c-di-GMP at a rate k3
% Clp.c-di-GMP > Clp + c-di-GMP at a rate k4
% Clp + PmanA > Clp.PmanA at a rate k5
% Clp.PmanA > Clp + PmanA at a rate k6
% Clp.PmanA > GFP at a rate K
% additionally, Clp is made at a rate k7 and decayed at a rate k8
% and c-di-GMP is made at a rate k9 and decayed at a rate k10

% Inputs: 
%      t: time variable
%      x: state variables as a vector
%     mu: concentration of DSF
%    clp: production rate of Clp
%
% Outputs:
%      dx_dt: RHS of the system of ODEs
%
% Usage:
% [t, x] = ode45(@final_setup_dsf,[time interval],[initial conditions],[],mu,clp);
% [t, x] = ode45(@final_setup_dsf,[0,1200],[4.98,0,2,2,0,0,0],[],mu,Clp);

%Parameters:
k1 = 1/60;        % formation of RpfG[P]                
k2 = 1/60;        % breaking down of c-di-GMP by RpfG[P]          
k3 = 2/60;        % binding rate of Clp + c-di-GMP           
k4 = 7/60;        % dissociation rate of Clp.c-di-GMP                
k5 = 5/60;        % binding rate of Clp and PmanA
k6 = 0.07/60;     % dissociation rate of Clp.PmanA
k7 = clp;         % production rate of Clp
k8 = 0.001/60;    % degradation rate of Clp
k9 = 0.01/60;     % production rate of c-di-GMP
k10 = 0.001/60;   % degradation rate of c-di-GMP
K = 1/60;         % maximal rate of GFP expression per promoter

D = mu;           % DSF (signal)
P = 0.083;        % PmanA (promoter) 

% vector of variables
% x(1) = RpfG
% x(2) = RpfG[P]
% x(3) = c-di-GMP
% x(4) = Clp
% x(5) = Clp.c-di-GMP
% x(6) = Clp.PmanA
% x(7) = GFP

%ODE system
%dx_dt(1) = -k1*D*x(1) + k2*x(3)*x(2) is the ODE for RpfG etc.

dx_dt(1) = -k1*D*x(1) + k2*x(3)*x(2);
dx_dt(2) = k1*D*x(1) - k2*x(3)*x(2);
dx_dt(3) = -k2*x(3)*x(2) - k3*x(3)*x(4) + k4*x(5) + k9 - k10*x(3);
dx_dt(4) = -k3*x(3)*x(4) + k4*x(5) - k5*P*x(4) + k6*x(6) + k7 - k8*x(4);
dx_dt(5) = k3*x(3)*x(4) - k4*x(5);
dx_dt(6) = k5*P*x(4) - k6*x(6);
dx_dt(7) = K*x(6);

%dx_dt is a column vector of ODEs

dx_dt=dx_dt';

%we transpose to allow MATLAB to solve the ODE system as a row vector

end

