clear all close all clc %% parameters alpha=1.5; %birth rate of uninfected population C=100; %2*10^6; %carrying capacity %k=6.7*10^(-8); lambda=3; %infected population death rate mu=20; %virus death rate q=10^(-6); %probability B=14.9254*10; %% phenotype discretization a=0; b=9; h=0.01; n=(b-a)/h; x=[a:h:b]; %% k as linear function aa=0.02; %0.002; k=aa*x; R0=(B-1)./(mu./(k*C)); pos=find(x==1); %% time discretization t0=0; tf=1000; dt=0.01; t=[0:dt:tf]; tmax=length(t); %% initial conditions p=zeros(n+1,1); i=zeros(n+1,1); %virus initial condition p(pos-5:pos+5) = 0.3; %infected cells initial condition i(pos-5:pos+5) = 0.3; %susceptible cells initial condition s = C; %% Boundary conditions p(1)=p(2); p(n+1)=p(n); i(1)=q*(i(2)-i(1))/h; i(n+1)=q*(i(n+1)-i(n))/h; %infected cells initial distribution: plot figure(1) plot(x, i, 'k', 'LineWidth',1.5) title('Initial distribution of infected cells','fontsize', 15) xlabel('r','fontsize',15) ylabel('cells','fontsize',15) %% method % inizialization tp=t0; iter=0; i_new=i; p_new=p; % save solutions i_pop=i; s_pop=s; i_int=sum(i_pop); for time=1:tmax% Time Loop iter=iter+1; for l=2:n; % Space Loop p_new(l)=p(l)+dt*(-k(l)*p(l)*s-mu*p(l)+B*lambda*i(l)); i_new(l)=i(l)+dt*(k(l)*p(l)*s-lambda*i(l))+(q*dt/(h^2))*(i(l+1)-2*i(l)+i(l-1)); end p_new(1)=p_new(2); p_new(n+1)=p_new(n); i_new(1)=q*(i_new(2)-i_new(1))/h; i_new(n+1)=q*(i_new(n+1)-i_new(n))/h; i=i_new; p=p_new; s=s+dt*(alpha*s*(1-(s/C)-(1/C)*(sum(i)))-sum(s*k'.*p)); int=sum(i); if mod (time,50)==0 i_pop=[i_pop, i]; s_pop=[s_pop, s]; i_int=[i_int, int]; tp=[tp t(time)]; end end figure(2) mesh(x,tp,i_pop') colorbar axis([0 9 0 1000]) %title('Distribution of infected cells in the viral phenotype space','fontsize', 15) xlabel('r','fontsize',15) ylabel('Time','fontsize',15) figure(3) plot(tp, s_pop, 'k', 'LineWidth',1.5) title('Dynamics of the susceptible population','fontsize', 15) xlabel('Time','fontsize',15) ylabel('Concentration of cells','fontsize',15) figure(4) plot(x, i_pop(:,end), 'k', 'LineWidth',1.5) title('Distribution of infected cells at t=1000','fontsize', 15) xlabel('r','fontsize',15) ylabel('Concentration of bacteria','fontsize',15) figure(5) plot(tp, i_int, 'k', 'LineWidth', 1.5) title('Dynamics of the total infected population','fontsize', 15) xlabel('Time','fontsize',15) ylabel('Concentration of cells','fontsize',15)