************************************************************
(1) c_agents_immune_oned.m
************************************************************
theta_1=2*Du*tau/delta^2;
theta_2=2*Di*tau/delta^2;

theta_T=2*DT*tau/delta^2;
nu=2*chiT*tau*cstar/delta^2;
Dc_ag=Dc*tau/delta^2;

T=T_pde/tau;
griglia=L/delta; %numeri pari
griglia_uninf=cell_uninf/delta; %numero punti di griglia con cellule non infette
griglia_inf=cell_inf/delta;     %numero punti di griglia con cellule infette
n_uninf=h_uninf*delta;
n_inf=h_inf*delta;

U=zeros(2*griglia+1,T_pde+1);
I=zeros(2*griglia+1,T_pde+1);
Tcell=zeros(2*griglia+1,T_pde+1);
c=zeros(2*griglia+1,T_pde+1);

vessel=ones(2*griglia+1,1);

% Maschera inflow: solo 40% sinistro del dominio
mask = zeros(2*griglia+1,1);
mask(1:round(0.4*length(mask))) = 1;

U(:,1)=[zeros(griglia-griglia_uninf,1)
      n_uninf*ones(2*griglia_uninf+1,1) 
      zeros((griglia-griglia_uninf),1)];
I(:,1)=[zeros(griglia-griglia_inf,1)
      n_inf*ones(2*griglia_inf+1,1) 
      zeros((griglia-griglia_inf),1)];

U_temp=U(:,1);
I_temp=I(:,1);
Tc_temp=Tcell(:,1);
c_temp=c(:,1);

for i=2:T+1
   % pressione = densita'  tumorale totale
   P = min(K, U_temp + I_temp);

   % --- MIGRAZIONE TUMORALE A PRESSIONE (sostituisce smorzata) ---
   U_temp = passeggiata_unod_pressione(U_temp, P, theta_1, K);
   I_temp = passeggiata_unod_pressione(I_temp, P, theta_2, K);

   % --- MOVIMENTO IMMUNE (chemotassi + ostacolo) ---
   Tc_temp = passeggiata_unod_chemotaxisdiffusion_obstacle(Tc_temp, c_temp, P, theta_T, nu, cstar, w_max);

   % --- NASCITA/MORTE ---
   U_temp = nascitamorte_pressione(U_temp, P, K*delta, p*tau, 0);
   I_temp = nascitamorte_pressione(I_temp, P, K*delta, 0, q*tau);
   Tc_temp= nascitamorte_pressione(Tc_temp,P,K*delta,0,qT*tau);

   % --- INFLOW IMMUNE condizionato alla maschera ---
   Tc_temp = inflow_masked(Tc_temp, alphaz * sum(I_temp) * delta * tau, vessel, mask);

   % --- INFEZIONE ---
   [U_temp,I_temp] = infezione_dued(U_temp, I_temp, beta*tau/delta);

   % --- KILLING IMMUNE ---
   U_temp = T_kill(U_temp, Tc_temp, min(1,zU*tau/delta));
   I_temp = T_kill(I_temp, Tc_temp, min(1,zI*tau/delta));

   % --- CHEMOATTRATTANTE ---
   c_temp = dinamicac_unod(c_temp,Dc_ag,I_temp,alphac*tau/delta,U_temp,gammac*tau/delta,cstar,qc*tau);

   if any(isnan(U_temp)) || any(isnan(I_temp)) || any(isnan(Tc_temp))
        i
   end
    
   % salvataggi
   t=(i-1)*tau;
   if t==floor(t)
        U(:,t+1)=U_temp;
        I(:,t+1)=I_temp;
        Tcell(:,t+1)=Tc_temp;
        c(:,t+1)=c_temp;
        
        figure(1)
        clf
        subplot(1,2,1)
        plot(U_temp/delta)
        hold on
        plot(I_temp/delta)
        subplot(1,2,2)
        plot(Tc_temp/delta)
   end
end




****************************************************************************************************
(2) c_main_immune_oned.m
****************************************************************************************************
clearvars
c_parameters_immune_oned

L=10;
T_pde=1000;
maxit=1;
delta=0.05; %mm
tau=0.01; %h
dx=0.05;
dt=1e-3;
cell_uninf=Ru; %larghezza iniziale dominio non infetto
cell_inf=Ri;   %largezza iniziale dominio infetto
h_uninf=0.9*K;
h_inf=0.1*K;
%save('c_parameters')

space=[-L L];

chiT=chiT/10;
zU=zU/2;
zI=zI/2;
q=q/5;
S_z=0.05;
w_max = 0.4 * K;


tic
c_pde_immune_oned
time_pde=toc
% save('pde','u1','u2','u3','u4','x','t','dx','dt','time_pde')
%%
tic
seed=1;
rng(seed)
c_agents_immune_oned
time_agents=toc
% save('agents'+string(seed),'U','I','Tcell','c','seed','time_agents')

%%

U_plot=U;
I_plot=I;
Tcell_plot=Tcell;
c_plot=c;

h = figure(1); clf
vid=VideoWriter('video');
% vid=VideoWriter('/fred/oz055/david/2023/immune/v_${SLURM_JOB_NAME}_${SLURM_ARRAY_JOB_ID}_${SLURM_ARRAY_TASK_ID}');
vid.FrameRate=4;
c_video_oned_immune


figure(2)
clf
U_sum=sum(U);
I_sum=sum(I);
T_sum=sum(Tcell);

u1_sum=sum(u1')*dx;
u2_sum=sum(u2')*dx;
u3_sum=sum(u3')*dx;

plot(0:T_pde,U_sum,'LineWidth',2)
hold on
plot(0:T_pde,I_sum,'LineWidth',2)
hold on
plot(0:T_pde,U_sum+I_sum,'LineWidth',2)
hold on
plot(0:T_pde,T_sum,'LineWidth',2)
hold on
plot(0:T_pde,u1_sum,'k:','LineWidth',2)
hold on
plot(0:T_pde,u2_sum,'k:','LineWidth',2)
hold on
plot(0:T_pde,u1_sum+u2_sum,'k:','LineWidth',2)
hold on
plot(0:T_pde,u3_sum,'k:','LineWidth',2)
axis([0 T_pde 0 max(u1_sum+u2_sum)])

legend('Uninfected','Infected','Total','Immune','Location','northwest')




***************************************************************************************************************
(3) c_parameters_immune_oned.m
***************************************************************************************************************
p=log(2)/37;
q=1/24;
alpha=3500; %virus rilasciati da una cellula
beta=7*1e-10*alpha/24; %beta Friedmanetal06 * alpha (vedi paper per il motivo)
Du=((6.2-2.6)/2/40/24)^2/p; %ogni unita' spaziale vale 1 mm
Di=Du;

K=1e4;
Ru=2.6; %raggio iniziale non infette in mm
Ri=1; %raggio iniziale infette in mm

beta=beta*1e6/K;

Dc=0.8/24; %MCK04
qc=2/24; %Emma
% qc=log(2)/0.5; %too much
cstar=5^(1/3); %mug/mm^2
alphac=2.5/K; %valore massimo di Adrianne diviso...; 1/h/cell
gammac=alphac/50;

qT=0.18/24; %OK
alphaz=qT/50; %arbitrario

zU=10/K/10;
zI=zU;

%DT=0.025/24; %Ats+20
DT=0.0023; %from microscopic considerations
chiT=1.65;
w_max = 0.4 * K;  % soglia di densita'  tumorale oltre cui la mobilita'  T si annulla




*****************************************************************************************************************
(4) c_pde_immune_oned.m
*****************************************************************************************************************
T=T_pde/dt;
x=-L:dx:L;
Nx=2*L/dx;

Nx_u=2*Ru/dx;
Nx_i=2*Ri/dx;

u1=zeros(Nx+1,T_pde+1);
u2=u1; u3=u1; u4=u1;

u=[zeros((Nx-Nx_u)/2,1)
      h_uninf*ones(Nx_u+1,1) 
      zeros((Nx-Nx_u)/2,1)];
i=[zeros((Nx-Nx_i)/2,1)
      h_inf*ones(Nx_i+1,1) 
      zeros((Nx-Nx_i)/2,1)];
z=zeros(Nx+1,1);
c=zeros(Nx+1,1);

u1(:,1)=u; u2(:,1)=i; u3(:,1)=z; u4(:,1)=c;

for k=1:T 
    rho = (u+i);

    % --- DIFFUSIONE TUMORALE: schema pressure-driven upwind ---
    u = update_tumor_pde_pressure(u, Du, K, p, dt, dx);
    i = update_tumor_pde_pressure(i, Di, K, 0, dt, dx);  % no crescita per infette

    % --- DIFFUSIONE CHEMOATTRATTANTE ---
    vicino_destra=[c(2:Nx+1);0];
    vicino_sinistra=[0;c(1:Nx)];
    c=c+Dc*(vicino_sinistra+vicino_destra-2*c)*dt/(dx^2);
    c(1)=c(2); c(Nx+1)=c(Nx);

    % --- MOVIMENTO IMMUNE: chemotassi upwind + diffusione ---
    c_destra=[c(2:Nx+1);0];
    z_destra=[z(2:Nx+1);0];
    w=chiT*(c_destra-c)/dx;
    Fz=max(w,0).*z - max(-w,0).*z_destra - DT*(z_destra-z)/dx;
    z=z+(dt/dx)*([0;Fz(1:Nx)]-[Fz(1:Nx);0]);

    % --- REAZIONI ---
    u=u+dt*( -beta*u.*i - zU*u.*z );
    i=i+dt*( beta*u.*i - q*i - zI*i.*z );
    z=z+dt*(-qT*z+alphaz*sum(i)*dx);
    c=c+dt*((alphac*i+gammac*u).*(cstar-c)-qc*c);

    if (mod(k,1/dt)==0)
        u1(:,k*dt+1)=u;
        u2(:,k*dt+1)=i;
        u3(:,k*dt+1)=z;
        u4(:,k*dt+1)=c;
        
        figure(1)
        clf
        sgtitle('t='+string(k*dt)+' h')
        subplot(1,2,1)
        plot(x,u,'LineWidth',1)
        hold on    
        plot(x,i,'LineWidth',1)
        hold on
        plot(x,c*200,'LineWidth',1)
        legend('u','i')
        
        subplot(1,2,2)
        plot(x,z,'LineWidth',1)
        hold on
        plot(x,c*10,'LineWidth',1)
        legend('z','c')
    end
end

t=0:T_pde;
u1=u1'; u2=u2'; u3=u3'; u4=u4';



***********************************************************************************
(5) c_video.m
***********************************************************************************
open(vid)

theta_1=2*Du*tau/delta;
theta_2=2*Di*tau/delta;
minAx=0;
maxAx=K;
griglia=L/delta; %numeri pari

vu=sqrt(Du*p/2);
vi=sqrt(Di*(beta*K-q)/2);

for i=1:10:T_pde+1
    clf
    plot((0:griglia)*delta,U_plot(:,i)/delta,'LineWidth',1)
    hold on
    plot(x,u1(i,:),'--')
    hold on
    yline(q/beta)
    hold on
    xline(L/2+Ru+(i-1)*vu,'--')
    hold on
    xline(L/2-Ru-(i-1)*vu,'--')
    hold on
    plot((0:griglia)*delta,(U_plot(:,i)+I_plot(:,i))/delta)
    plot((0:griglia)*delta,I_plot(:,i)/delta,'LineWidth',1)
    hold on
    plot(x,u2(i,:),'--')
    hold on
    yline(p*(beta*K-q)/beta/(beta*K+p))
    axis([0 griglia*delta minAx maxAx])
    title('Cellule infette', 'q='+string(q)+', \beta='+string(beta)+', \theta_2='+string(theta_2)+...
        ', \delta='+string(delta)+', \tau='+string(tau))
    legend('Agenti','PDE')

    frame=getframe(h);
    writeVideo(vid,frame);
end


close(vid)




********************************************************************************************
(6) c_video_oned_immune.m
********************************************************************************************
open(vid)

minAx=0;
maxAx=K;

for i=1:10:T_pde+1
    clf
    sgtitle('t='+string(i-1)+' h')
    subplot(1,2,1)
    plot(-L:delta:L,U_plot(:,i)/delta,'LineWidth',1)
    hold on    
    plot(-L:delta:L,I_plot(:,i)/delta,'LineWidth',1)
    hold on
    plot(x,u1(i,:),'k:','LineWidth',1)
     hold on
    plot(x,u2(i,:),'k:','LineWidth',1)

    subplot(1,2,2)
    plot(-L:delta:L,Tcell_plot(:,i)/delta,'LineWidth',1)
    hold on
    plot(-L:delta:L,c_plot(:,i)/10,'LineWidth',1)
    hold on
    plot(x,u3(i,:),'k:','LineWidth',1)
    hold on
    plot(x,u4(i,:)/10,'k:','LineWidth',1)

    plot(-L:delta:L,c_plot(:,i),'LineWidth',1)
    hold on
    plot(x,u4(i,:),'k:','LineWidth',1)

    frame=getframe(h);
    writeVideo(vid,frame);
end


close(vid)






***********************************************************************************
(7) dinamicac_unod.m
***********************************************************************************
function c=dinamicac_unod(c,D,I,alphac,U,gammac,cstar,qc)
%Funzione che prende in input un vettore colonna a valori interi positivi 
%di lunghezza almeno 3 e restituisce il vettore ottenuto a seguito di
%diffusione, decadimento e produzione della popolazione I. 
%Condizioni al bordo di Neumann omogenee

c_right=[c(2:end);0];
c_left=[0;c(1:end-1)];

c=c+D*(c_left+c_right-2*c)+(alphac*I+gammac*U).*(cstar-c)-qc*c;
c(1)=c(2);
c(end)=c(end-1);





************************************************************************************
(8) infezione_dued.m
************************************************************************************
function [U,I] = infezione_dued(U, I, beta)
%Funzione che prende in input due matrici U, I a valori interi positivi, una
%probabilità di infezione beta e restituisce la matrice ottenuta dopo
%l'infezione


prob_max=1;

A=binornd(U,min(prob_max, beta*I));
U=U-A;
I=I+A;





**************************************************************************************
(9) inflow_masked.m
**************************************************************************************
function v = inflow_masked(v, p_inflow, vessel, mask)
% Inflow condizionato a regioni definite da 'mask'
% mask vale 1 dove l'inflow e' permesso, 0 altrove
% vessel vale 1 nei punti dei vasi

a = v(:);
b = vessel(:) .* mask(:);  % solo dove entrambe sono 1

a(b == 1) = a(b == 1) + (rand(sum(b), 1) < p_inflow);

v(:) = a;




**************************************************************************************
(10) nascitamorte_pressione.m
**************************************************************************************
function v = nascitamorte_pressione(v, P, K, p, q)
%Funzione che prende in input una matrice v a valori interi positivi, una
%matrice pressione P delle stesse dimensioni, una carrying capacity K, una 
%probabilità di nascita p e una di morte q e restituisce la matrice 
%ottenuta dopo riproduzione e morte. Se la pressione è superiore a K, le
%cellule muoiono prima a tasso p per riequilibrare.
%NB: le nascite precedono le morti, muoiono anche i neonati


v=v+binornd(v,p*(1-P/K).*(P<K));
v=v-binornd(v,p*(P/K-1).*(P>=K)); %muoiono le eventuali cellule in eccesso
v=binornd(v,1-q);




****************************************************************************************
(11) passeggiata_unod_chemotaxisdiffusion_obstacle.m
****************************************************************************************
function w = passeggiata_unod_chemotaxisdiffusion_obstacle(v, c, P, theta, chi, cmax, w_max)
% Movimento stocastico con chemotassi e ostacolo all'infiltrazione
% v     = T-cell attuali
% c     = chemoattrattante
% P     = densita'  tumorale infetto + non infetto (ostacolo)
% theta = probabilita'  massima di diffusione
% chi   = coefficiente di chemotassi
% cmax  = scala normalizzante per chemotassi
% w_max = soglia oltre la quale si blocca il movimento

n = length(v);
w = zeros(n,1);

% Calcolo della mobilita'  locale
densita = P;
mobilita = max(0, 1 - densita / w_max);

% Movimento interno (escludendo i bordi)
theta_eff = theta * mobilita(2:end-1);
chi_eff = chi * mobilita(2:end-1);

prob_s = theta_eff / 2 + chi_eff .* subplus(c(1:end-2) - c(2:end-1)) / (2 * cmax);
prob_d = theta_eff / 2 + chi_eff .* subplus(c(3:end) - c(2:end-1)) / (2 * cmax);
prob_c = 1 - prob_s - prob_d;

% Multinomial sampling per celle interne
b = mnrnd(v(2:end-1), [prob_s, prob_c, prob_d]);

% Ricomposizione del nuovo vettore
w = [b(:,1); 0; 0] + [0; b(:,2); 0] + [0; 0; b(:,3)];

% BORDO SINISTRO
theta_edge = theta * mobilita(1);
chi_edge = chi * mobilita(1);
p_dx = theta_edge / 2 + chi_edge * subplus(c(2) - c(1)) / (2 * cmax);
p_cx = 1 - p_dx;

a = mnrnd(v(1), [p_cx, p_dx]);
w(1) = a(1);
w(2) = w(2) + a(2);

% BORDO DESTRO
theta_edge = theta * mobilita(end);
chi_edge = chi * mobilita(end);
p_sx = theta_edge / 2 + chi_edge * subplus(c(end-1) - c(end)) / (2 * cmax);
p_dx = 1 - p_sx;

a = mnrnd(v(end), [p_sx, p_dx]);
w(end-1) = w(end-1) + a(1);
w(end) = a(2);




****************************************************************************************
(12) passeggiata_unod_pressione.m
****************************************************************************************
function w = passeggiata_unod_pressione(v,P,theta,K)
% Passeggiata unidimensionale guidata da pressione
% v = vettore popolazione
% P = pressione locale (stesso size)
% theta = probabilita'  massima di moto
% K = carrying capacity

n=length(v);
P=min(P,K);

prob_s = theta*subplus(P(2:n-1)-P(1:n-2))/(2*K);
prob_d = theta*subplus(P(2:n-1)-P(3:n))/(2*K);

b = mnrnd(v(2:n-1), [prob_s 1-prob_s-prob_d prob_d]);

w = [b(:,1);0;0] + [0;b(:,2);0] + [0;0;b(:,3)];

% Estremi riflettenti
a = binornd(v(1), theta*subplus(P(1)-P(2))/(2*K));
w(1)=w(1)+v(1)-a; w(2)=w(2)+a;

a = binornd(v(n), theta*subplus(P(n)-P(n-1))/(2*K));
w(n)=w(n)+v(n)-a; w(n-1)=w(n-1)+a;
end



****************************************************************************************
(13) T_kill.m
****************************************************************************************
function v = T_kill(v, T_cell, z)
%Funzione che prende in input una matrice v a valori interi positivi, una
%matrice T_cell delle stesse dimensioni, una probabilità z di uccisione da 
%parte del sistema immunitario e  restituisce la matrice 
%ottenuta dopo l'azione del sistema immunitario

prob_max=1;

v=binornd(v,1-min(prob_max,z*T_cell));





****************************************************************************************
(14) update_tumor_pde_pressure.m
****************************************************************************************
function u = update_tumor_pde_pressure(u, D, K, p, dt, dx)
%UPDATE_TUMOR_PDE_PRESSURE
% Aggiorna la densita'  tumorale u secondo diffusione con pressione
% (upwind scheme) e crescita logistica.
%
% u  : densita'  tumorale (Nx+1 x 1)
% D  : coefficiente di diffusione
% K  : capacita'  portante
% p  : tasso di proliferazione
% dt : passo temporale
% dx : passo spaziale

Nx = length(u) - 1;

% Calcolo flussi (upwind per la pressione)
u_destra   = [u(2:end); 0];
u_sinistra = [0; u(1:end-1)];

% Flusso con upwind: diffusione modulata dalla densita' 
Fu = -D * (u_destra - u) / dx;

% Aggiornamento tramite divergenza dei flussi
u = u + (dt/dx) * ([0; Fu(1:Nx)] - [Fu(1:Nx); 0]);

% Reazione: crescita logistica sotto pressione
rho = u;
u = u + dt * (p * (1 - rho./K) .* u);

% Condizioni al bordo di Neumann (no-flux)
u(1)     = u(2);
u(end)   = u(end-1);

end
