function [ dxdt ] = HumanBodyFB(t,x)
%%%%%%%%%%%%%%%%%%%%Variable initialization%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Flow rate blood %%%     L/min
Qc = 5.58; %L/min
Qlung = Qc;
Qvenous = Qc;
Qarterial = Qc;
Qliver = 0.25 * Qc;
Qliverin = Qliver * 0.25;
Qintestine = Qliver * 0.75;
QintestineHealthy = Qintestine * 0.90;
QintestineSick = Qintestine * 0.10;
Qkidney = 0.19 * Qc;
Qrpt = 0.18 * Qc;
Qspt = 0.38 * Qc;

%%% Tissue Volumes %%%      L
BV = 70;
Vlung = 0.008 * BV;
Vvenous = 0.0557 * BV;
Varterial = 0.0243 * BV;
Vliver = 0.0260 * BV;
Vkidney = 0.0044 * BV;
Vrpt = 0.0516 * BV;
Vspt = 0.76 * BV;
Vintestine = 0.07 * BV;
VintestineHealthy = 0.9 * Vintestine;
VintestineSick = 0.1 * Vintestine;

%Amount of cells in human body
%HBcells = 10^14;
%Avogadro = 6,02214 * 10^23;
%Example:
%RTrpt = ((100000 * HBcells) / BV) / Avogadro;
%%% Receptor total %%%      uM
RTrpt = 0.237;
RTspt = 0.237;
RTliver = 0.237;
RTkidney = 0.237;
RTvenous = 0.237;
RTlung = 0.237;
RTarterial = 0.237;
RTintestineHealthy = 0.237;
RTintestineSick = 10.716;

%%% VLP Elimination rate %%%    uM
Kel = 0.000143;					%Biodistribution studies of protein cage nanoparticles demonstrate broad tissue distribution and rapid clearance in vivo

%%% VLP Renal removal rate %%%  uM
Kkidney = 0.0000403;				%Biodistribution studies of protein cage nanoparticles demonstrate broad tissue distribution and rapid clearance in vivo

%%% med Elimination rate %%%    uM
Kelmed = 0.05;					%Calculated by half-life from wiki...

%%% med Renal removal rate %%%  uM
%Kkidneymed = Kelmed * 1000000;

%%% packeging constant %%%      -
ymedvlp = 300;					%Calculated by mass

%%% Affinity constant for receptors %%%		http://www.bindingdb.org/bind/index.js
K = 0.001; %Nog keer 1 miljoen

%%% Uptake for MED to organs %%%
Pblungmed = 100;
Pblivermed = 100;
Pbkidneymed = 100;
Pbrptmed = 100;
Pbsptmed = 100;
PbintestineHealthymed = 100;
PbintestineSickmed = 100;
Pbvenous = 100;
Pbarterial = 100;

%%%%%%%%%%%%%%%%%%%%Receptor FREE concentration%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RFrpt = RTrpt - (x(1) - x(21));
RFspt = RTspt - (x(2) - x(22));
RFliver = RTliver - (x(3) - x(23));
RFkidney = RTkidney - (x(4) - x(24));
RFvenous = RTvenous - (x(5) - x(25));
RFlung = RTlung - (x(6) - x(26));
RFarterial = RTarterial - (x(7) - x(27));
RFintestineHealthy = RTintestineHealthy - (x(8) - x(28));
RFintestineSick = RTintestineSick - (x(9) - x(29));

%%%%%%%%%%%%%%%%%%%%Explanation of the dxdt%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dxdt(1) = Crpt
% dxdt(2) = Cspt
% dxdt(3) = Cliver
% dxdt(4) = Ckidney
% dxdt(5) = Cvenous
% dxdt(6) = Clung
% dxdt(7) = Carterial
% dxdt(8) = IntestineHealthy
% dxdt(9) = IntestineSick

%Total VLP Concentration
dxdt(1) = Qrpt / Vrpt * ( x(27) - x(21)) - Kel * x(1);
dxdt(2) = Qspt / Vspt * ( x(27) - x(22)) - Kel * x(2);
dxdt(3) = Qliverin / Vliver * x(27) + QintestineHealthy / Vliver * x(28) - Qliver / Vliver * x(23) + QintestineSick / Vliver * x(29) - Kel * x(3);
dxdt(4) = Qkidney/Vkidney * (x(27) - x(24)) - Kel * x(4) - ( Kkidney * x(4) / Vkidney);

%Venous flow from all organs
Cj = [x(21), x(22), x(23), x(24)];
Qj = [Qrpt, Qspt, Qliver, Qkidney];

som=0;
for i = 1:length(Qj)
    som = som + Qj(i) * Cj(i);
end

dxdt(5) = 1 / Vvenous * ( som - Qvenous * x(25))- Kel * x(5);
dxdt(6) = Qlung / Vlung * ( x(25) - (x(26)) - Kel * x(6));
dxdt(7) = 1 / Varterial * ( Qlung * x(26) - Qarterial * x(27)) - Kel * x(7);
dxdt(8) = Qintestine  / VintestineHealthy * x(27) - QintestineHealthy  / VintestineHealthy * x(28) - QintestineSick / VintestineHealthy * x(28) - Kel * x(8);
dxdt(9) = QintestineSick / VintestineSick * x(28) - QintestineSick / VintestineSick * x(29) - Kel * x(9);

%%%%%%%%%%%%%%%%%%%%%%FREE VLPs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Total FREE VLP Concentration
dxdt(21) = Qrpt / Vrpt * ( x(27) - x(21)) - Kel * x(21) - K * x(21) * RFrpt;
dxdt(22) = Qspt / Vspt * ( x(27) - x(22)) - Kel * x(22) - K * x(22) * RFspt;
dxdt(23) = Qliverin / Vliver * x(27) + QintestineHealthy / Vliver * x(28) - Qliver / Vliver * x(23) + QintestineSick / Vliver * x(29) - Kel * x(23) - K * x(23) * RFliver;
dxdt(24) = Qkidney/Vkidney * (x(27) - x(24)) - Kel * x(24) - K * x(24) * RFkidney;

%Venous flow from all organs
Cj = [x(21), x(22), x(23), x(24)];
Qj = [Qrpt, Qspt, Qliver, Qkidney];

som=0;
for i = 1:length(Qj)
    som = som + Qj(i) * Cj(i);
end

dxdt(25) = 1 / Vvenous * ( som - Qvenous * x(25))- Kel * x(25) - K * x(25) * RFvenous;

dxdt(26) = Qlung / Vlung * ( x(25) - (x(26)) - Kel * x(26)) - K * x(26) * RFlung;
dxdt(27) = 1 / Varterial * ( Qlung * x(26) - Qarterial * x(27)) - Kel * x(27) - K * x(27) * RFarterial;
dxdt(28) = Qintestine  / VintestineHealthy * x(27) - QintestineHealthy  / VintestineHealthy * x(28) - QintestineSick / VintestineHealthy * x(28) - Kel * x(28) - K * x(28) * RFintestineHealthy;
dxdt(29) = QintestineSick / VintestineSick * x(28) - QintestineSick / VintestineSick * x(29) - Kel * x(29) - K * x(21) * RFrpt - K * x(29) * RFintestineSick;

%%%%%%%%%%%%%%%%%%%%%%MED PACK%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Organs according to standard flow calculations
%RPT blood flow
dxdt(11) = Qrpt / Vrpt * ( x(17) / Pbarterial - (x(11) / Pbrptmed )) + Kel * x(1) * ymedvlp;
%SPT blood flow
dxdt(12) = Qspt / Vspt * ( x(17) / Pbarterial - (x(12) / Pbsptmed )) + Kel * x(2) * ymedvlp;
%Healthy intestines from Fres blood flow
dxdt(18) = Qintestine  / VintestineHealthy * x(17) / Pbarterial - QintestineHealthy  / VintestineHealthy * x(18) / PbintestineHealthymed - QintestineSick / VintestineHealthy * x(18) / PbintestineHealthymed + Kel * x(8) * ymedvlp;
%Tumor intestine from Healthy intestines
dxdt(19) = QintestineSick / VintestineSick * x(18) / PbintestineHealthymed - QintestineSick / VintestineSick * x(19) / PbintestineSickmed + Kel * x(9) * ymedvlp;
%Liver blood flow from Healthy intestine and from Fresh blood flow
dxdt(13) = Qliverin / Vliver * x(17) / Pbarterial + QintestineHealthy / Vliver * x(18) / PbintestineHealthymed - Qliver / Vliver * x(13) / Pblivermed + QintestineSick / Vliver * x(19) / PbintestineSickmed + Kel * x(3) * ymedvlp;
%Kidney from Healthy blood flow, also removes some of the medicine
dxdt(14) = Qkidney/Vkidney * (x(17) / Pbarterial - x(14) / Pbkidneymed ) - ( Kelmed * x(14) / Vkidney * Pbkidneymed ) + Kel * x(4) * ymedvlp;

%Venous flow from all organs
Cj = [x(11), x(12), x(13), x(14)];
Qj = [Qrpt, Qspt, Qliver, Qkidney];
Pbj = [Pbrptmed, Pbsptmed, Pblivermed, Pbkidneymed];

som=0;
for i = 1:length(Qj)
    som = som + Qj(i) * Cj(i) / Pbj(i);
end

dxdt(15) = 1 / Vvenous * ( som - (Qvenous * x(15) / Pbvenous)) + Kel * x(5) * ymedvlp;

%Lung from Venous
dxdt(16) = Qlung / Vlung * ( x(15) / Pbvenous - (x(16) / Pblungmed )) + Kel * x(6) * ymedvlp;

%Arterial from Lung
dxdt(17) = 1 / Varterial * ( Qlung * x(16) / Pblungmed - (Qarterial * x(17) / Pbarterial )) + Kel * x(7) * ymedvlp;

dxdt(10) = 0;
dxdt(20) = 0;

%transpose result
dxdt = dxdt';

end