% model is similar to model 12
% H(phi) in chemotaxis is changed to 0.5*(tanh(phi)+1)
t_max = 160; %90 150 500
x_1 = 0;
x_2 = 40;
delta_t = .01; %0.01 0.002
delta_x = 0.2; %0.2 0.1
wound_edge = 25;
wound_frac = wound_edge/(x_2-x_1);
healthy_frac = (x_2-wound_edge)/(x_2-x_1);
tot_space_steps = (x_2-x_1)/delta_x;
N0 = [repmat(0,wound_frac*tot_space_steps+1,1);repmat(1,healthy_frac*tot_space_steps,1)];
M0 = [repmat(1,wound_frac*tot_space_steps+1,1);repmat(0,healthy_frac*tot_space_steps,1)];
D0 = [repmat(0.1,wound_frac*tot_space_steps+1,1);repmat(1,healthy_frac*tot_space_steps,1)];
G0 = [repmat(0,wound_frac*tot_space_steps+1,1);repmat(0,healthy_frac*tot_space_steps,1)];
F0 = [repmat(0,wound_frac*tot_space_steps+1,1);repmat(1,healthy_frac*tot_space_steps,1)];
C0 = [repmat(0,wound_frac*tot_space_steps+1,1);repmat(1,healthy_frac*tot_space_steps,1)];
S0 = [repmat(1,wound_frac*tot_space_steps+1,1);repmat(1,healthy_frac*tot_space_steps,1)];
W0 = [repmat(0.5,wound_frac*tot_space_steps+1,1);repmat(0,healthy_frac*tot_space_steps,1)]; %initial bacterial innoculum
phi0 = [repmat(1,wound_frac*tot_space_steps+1,1);repmat(0,healthy_frac*tot_space_steps,1)];
H0 = [repmat(0,wound_frac*tot_space_steps+1,1);repmat(0,healthy_frac*tot_space_steps,1)];
A0 = [repmat(0,wound_frac*tot_space_steps+1,1);repmat(0,healthy_frac*tot_space_steps,1)];
N_bound = 1;
M_bound = 0;
D_bound = 1;
G_bound = 0;
F_bound = 1;
C_bound = 1;
S_bound = 1;
W_bound = 0;
phi_bound = 0;
H_bound = 0;
A_bound = 0;

chi_n = .001; %0.001
sigma_n = 1;
gamma_ncs = 1;
kappa_n = 1.2;
gamma_nc = 1;
delta_nc = 0; 
delta_mf = 0.2;
delta_ms = 0.02; 
delta_m = 0;
beta_gs = .16; 
delta_gc = 0;
beta_gf = .12; 
S_bar = 3; %3
beta_sc = .4; %0.4
D_n = 1;
D_f =.5; 
chi_f = .01; %0.1
sigma_f = 1;
rho_b = 0.6; 
beta_fb = .04; 
delta_sn = 0.4;
delta_sf = .4; 
beta_sm = 4; %6 or 4 for poor clotting
delta_s = 0.1; 
eta = 10; 
beta_fc = 0.8;
beta_ssf = 0.2;
kappa_f = .5; 
beta_fcs = 1; %3 or 1 cells in chronic wound are senescent so may struggle to uptake nutrients and grow
gamma_fcs = 1;
delta_fc = 0.4; %0.4
rho_f = 1;
D_c = 1.5;
beta_cb = 1; 
delta_cn = 0; 
delta_cf = 0.6; 
D_s = 1.5; %1.5
beta_sn = 0.24; 
beta_ssn = 0.2;
beta_sf = 0.8;
C_bar = 0.3; %.3 
gamma_gf = 0.1;

gamma_b = 0.85; %0.6 or 0.85 for poor vascuature

% bacteria params
alpha_dmw = 0.1; %.1
alpha_gmw = 0.1; %.1
alpha_fmw = .01; %.01 or 3 to make sure bacteria need to be killed for wound to heal
alpha_cmw = .4; %.4
gamma_gw = .5; %0.1 .5
beta_wc = .2; %.2
beta_wm = 0.2; %0.2
beta_wn = 0.001; %0.001
beta_sw = .4;
delta_w = 0.2;
delta_wf = 0.2; %0.2
beta_mn = 0.004; %0.004
beta_mg = .004; %0.004
beta_md = .04; %0.04
N_bar = 0.5;
delta_nw = 1; %0.04

%metal ion params
beta_nh = 0; 
beta_fh = 0; 
beta_sfh =10; % 10
D_h = 1.5; %1.5
delta_hf = .3; %.3
delta_hn = 0; 
delta_h = 0;

%stimulating granulation tissue production params
beta_gfh = 0; 


%metal ions vs bactera
D_a = 1.5; %1.5
delta_aw = .1; %.1
delta_a = .1; %.1
delta_wa = 2; %20

%dressing and wound depths
H_d=1;
H_w=1;

%ion release rate
delta_phi = 0.05; %0.05

%cell migration rates through fibres
sigma_nphi = 1; %1
sigma_fphi = 1; %1

%ion concentrations in fibres
alpha = 5; %known as gamma in paper, set to 0 for just antimicrobial ions

vol_frac=0.05;

n=1e-10;








%conditions for healed steady state in code
delta_nc = 1/(1+1/gamma_ncs)*(1-1/kappa_n)*(1+1/gamma_nc);
tmp = beta_fb*vasculture(1,0,rho_b,gamma_b,eta)-delta_fc/(1+(1+rho_f)/gamma_fcs);
kappa_f = beta_fc/(beta_fc+tmp);
delta_cn = beta_cb*vasculture(1,0,rho_b,gamma_b,eta)-delta_cf;
delta_s = beta_sn+beta_ssn+beta_sf+beta_ssf-delta_sn-delta_sf;
[N,M,D,G,F,C,S,W,phi,H,A] = solve_metal_ion_wound_healing_bacteria_system(t_max,x_1,x_2,delta_t,delta_x,N0,M0,D0,G0,F0,C0,S0,W0,phi0,H0,A0,N_bound,M_bound,D_bound,G_bound,F_bound,C_bound,S_bound,W_bound,phi_bound,H_bound,A_bound,chi_n,sigma_n,gamma_ncs,kappa_n,gamma_nc,delta_nc,delta_mf,delta_ms,delta_m,beta_gs,delta_gc,beta_gf,S_bar,beta_sc,D_n,D_f,chi_f,sigma_f,rho_b,beta_fb,delta_sn,delta_sf,beta_sm,delta_s,eta,beta_fc,beta_ssf,kappa_f,beta_fcs,gamma_fcs,delta_fc,rho_f,D_c,beta_cb,delta_cn,delta_cf,D_s,beta_sn,beta_ssn,beta_sf,C_bar,gamma_gf,alpha_gmw,alpha_dmw,alpha_fmw,alpha_cmw,beta_wc,beta_wm,beta_wn,beta_sw,delta_w,delta_wf,gamma_gw,N_bar, beta_mn,beta_mg,beta_md,delta_nw,beta_nh,beta_fh,beta_sfh,D_h,delta_hf,delta_hn,delta_h,D_a,delta_aw,delta_a,delta_wa,gamma_b,beta_gfh,sigma_nphi,sigma_fphi,delta_phi,alpha,H_d,H_w,vol_frac,n);

k = vasculture(D,G,rho_b,gamma_b,eta);
heav_fn = 0.5*(tanh(phi/n)+1);

no_of_iterations = t_max/8 + 1;
hold on
for i = 1:no_of_iterations
    solution = F(:,1+8*(i-1)/delta_t);
    plot(0:delta_x:x_2,solution)
end

%checks that the healed steady state conditions are satisfied
check_1 = beta_sn+beta_ssn+beta_sf+beta_ssf-delta_s-delta_sn-delta_sf;
check_2 = beta_cb/(gamma_b^eta+1)-delta_cn-delta_cf;
check_3 = 1/(1+1/gamma_ncs)*(1-1/kappa_n)-delta_nc/(1+1/gamma_nc);
check_4 = beta_fb/(gamma_b^eta+1)+beta_fc*(1-1/kappa_f)-delta_fc/(1+(1+rho_f)/gamma_fcs);