% Program tio2_10.m

% Global variables
global Dose r k13 k43 cl dr br k loadmax Ntotal phi rho nvar mcrit mcr2 lambda conc c dr2 k13_1 k13_2 alb ib sa satio2 Rec dr2 tfinal AMRec;
  clc;
% EXPOSURE INFORMATION-----------------------------------------------------------
% conc
% concentration (mg/m^3)
  conc = 10;
% starting time and end time for exposure period in days
 t0=0;
 tfinal=690;
% particle characteristics----------------------------------------------------
% density  (g/cm^3)
  density = 4.25;
% median aerodynamic diameter (um)
  diameter = 0.73;
  sa=60.67; % surface area (m^2/g)
  sa=sa*10^4/10^3; % (cm^2/mg)
%--------------------------------------------------------------------------
satio2=6.67; % surface area (m^2/g) of TiO2 as benchmark
satio2=satio2*10^4/10^3; % (cm^2/mg)
 
%**************************************************************************
% PHYSIOLOGICAL Parameters
% ventilation (m^3/day)
  ventilation= 0.15;        % (l/min)
  ventilation=ventilation/1000;    % conversion to m3/min
  ventilation=ventilation*60;      % conversion to m3/hour
  ventilation=ventilation*7;       % it is inhalation at 6 hours/day
  ventilation= ventilation; % (m^3/day)
% deposition fraction from MDDP programme
  depfract = 0.135;
% correction factor
  corfact = 0.714;
% Dose (mg/day)
D=conc*ventilation*depfract*corfact; % mg/day
%--------------------------------------------------------------------------
% Model parameters for EXPOSURE-DOSE Modelling
% (Phagocytosis rate day-1)
r = 1/0.25;
% (clearance day-1)
cl = 0.015;

% (interstitialized day-1)
% under normal circumstances
k13_1 = 0.03;
% during overload
k13_2 =1.8;
% (lymph day-1)
k43=0.01;
% Overload parameters
% rho is parameter from mobile to decay AMs (Stober)
rho=0.036;
% phi is rate of sequestration
phi=dr*1.0;


% lambda and c coefficients of for the retardation of clearance
lambda=5.8*(satio2/sa);
c=15;
  
% Model parameters for DOSE-RESPONSE modelling
% parameters for PMNs
  Rec=1.7/satio2;
  dr2=0.01; % death rate PMN
  
  % parameter for AMs
  AMtot=7.0;
  AMRec=0.7/sa;
% Death rate day-1, (Stober) for AM
  dr = 1.0/7;
%--------------------------------------------------------------------------
% Differential equations initial conditions

  Ntotal=11;
  tol=0.00001;
  trace=1;

  mcr2 = mcrit;
  N=1;
  Dose=D;
 
  tspan=[t0 tfinal];
  for i=1:(Ntotal-1)
     x0=0;
  end
  x0(Ntotal)=AMtot;
  
  options =odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6]);
  [time10,x]=ode45('model',tspan,x0,options);

% Post simulation conversion into appropriate units for plotting
% total burden
  burden10 = (x(:,1)+x(:,2)+x(:,3)+x(:,5)+x(:,6)+x(:,7)+x(:,8)+x(:,9)); % unit mg
  iburden=x(:,3)+x(:,7)+x(:,8)+x(:,9);
  lburden10=x(:,4);
  PMN=x(:,10);
  AM=x(:,11);