%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% BSc project       M.M. Mul                (mathias_mul@hotmail.com) %%%
%%% Coordinator       Dr.ir. F.J. Vermolen    (f.j.vermolen@tudelft.nl) %%%
%%% Period            february-june 2012                                %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Basic wound

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Model parameters %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%White blood cell parameters
white_blood_cell_size = 0.35;           % Diameter of white blood cell (space units) 
wbc_concentration = 0.5;                % Concentration of white blood cells in blood
active_range_wbc = 0.1;                 % Range in which virusses are annihilated (space units)
rvm_wbc = 0.1;                          % Random vertical movement of wbc inside blood vessel
beta = 1.5;                             % Mobility parameter of the white blood cell boundary
alpha = 5;                              % Contraction coefficient white blood cells
 
%Wound parameters
wound_center =[2,-6];                   % Center of wound
wound_horizontal_size =5;               % Horizontal size of wound (space units) 
wound_vertical_size =2;                 % Vertical size of wound (space units) 

%Blood vessel parameters
L_blood_vessel = 24;                    % The length of the blood vessel (integer space units)
thickness_blood_vessel = 2;             % Thickness of blood vessel (space units) 
thickness_blood_vessel_cell = 0.3;      % Thickness of blood vessel cell (space units) 
stiffness_bvc = 20;                     % Stiffness coefficient of blood vessel cells

%Virus colony parameters
virus_density = 0.8;                    % Virus density at t=0
sigma_virus = 0.8;                      % Mobility of virusses
D = 2;                                  % Diffusion coefficient of virusses 
br_virus = 0.8;                         % Birth rate of virus population
dr_virus = 30;                          % Death rate of virus population

%Heartbeat parameters
Blood_flow = 4.5;                       % Blood flow velocity (velocity units)
cc = 2;                                 % Length of cardiac cycle (time units)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Plot Window %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Complete situation
plot_window = [-L_blood_vessel/2,L_blood_vessel/2,min(-5,wound_center(2)-1.5*wound_vertical_size),max(5,wound_center(2)+1.5*wound_vertical_size)];
%Close-up on the wound
% plot_window = [wound_center(1)-wound_horizontal_size,wound_center(1)+wound_horizontal_size,wound_center(2)-wound_vertical_size,wound_center(2)+wound_vertical_size];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Time discretization %%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = 200;                        % calculations are made on time interval [0,T] 
dt = 0.1;                       % time interval
K = T/dt+1;                     % total number of timesteps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Blood vessel cells %%%%%%
%%%%%%%% & Blood vessel %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N_bvc = 12;                                                     % Number of nodes on blood vessel cell
L_blood_vessel_cell = 2;                                        % Length of blood vessel cell
N_blood_vessels = ceil(L_blood_vessel/L_blood_vessel_cell);     % Number of blood vessels on one side
%Defining starting position blood vessel cell centers
xc_bvc_up = (0:L_blood_vessel_cell:L_blood_vessel)-L_blood_vessel/2+L_blood_vessel_cell/2;
yc_bvc_up = ones(1,N_blood_vessels+(L_blood_vessel/L_blood_vessel_cell~=0))*(1/2)*(thickness_blood_vessel+thickness_blood_vessel_cell);
xc_bvc_down = xc_bvc_up;
yc_bvc_down = ones(1,N_blood_vessels+(L_blood_vessel/L_blood_vessel_cell~=0))*(-1/2)*(thickness_blood_vessel+thickness_blood_vessel_cell);
% Defining the shape of the blood vessel cells
x_initial_bvc = (1/2)*L_blood_vessel_cell*cos(((0:N_bvc-1)')*2*pi/N_bvc)';
y_initial_bvc = (1/2)*thickness_blood_vessel_cell*sin(((0:N_bvc-1)')*2*pi/N_bvc)'; 
%Defining starting postition nodes on cell boundary of blood vessel cells
x_bvc_up = ones(N_blood_vessels+(L_blood_vessel/L_blood_vessel_cell~=0),1)*x_initial_bvc + (ones(N_bvc,1)*xc_bvc_up)';   
y_bvc_up = ones(N_blood_vessels+(L_blood_vessel/L_blood_vessel_cell~=0),1)*y_initial_bvc + (ones(N_bvc,1)*yc_bvc_up)';
x_bvc_down = ones(N_blood_vessels+(L_blood_vessel/L_blood_vessel_cell~=0),1)*x_initial_bvc + (ones(N_bvc,1)*xc_bvc_down)';   
y_bvc_down = ones(N_blood_vessels+(L_blood_vessel/L_blood_vessel_cell~=0),1)*y_initial_bvc + (ones(N_bvc,1)*yc_bvc_down)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Heart beat %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_flow = Blood_flow;          %Blood flow velocity
min_flow = max_flow*(2/3);
lsf = (1/4);                    %length of the systolic phase of one heart beat [length systolic phase + length diastolic phase=1]
%Calculation of periodic heart beat function using splines for interpolation.
hbt_alpha = 0;
hbt_beta = cc*lsf;
hbt_gamma = cc;
x_hb_1 = 0:dt*lsf:hbt_beta;
x_hb_2 = hbt_beta:dt*lsf:cc;
t_hb_1=(1/(hbt_beta-hbt_alpha))*(x_hb_1-hbt_alpha);
t_hb_2=(1/(cc-hbt_beta))*(x_hb_2-hbt_beta);
f_hb_1 =(1-t_hb_1).*min_flow+t_hb_1.*max_flow+t_hb_1.*(1-t_hb_1).*(-(max_flow-min_flow).*(1-t_hb_1)+(max_flow-min_flow).*t_hb_1);
f_hb_2 =(1-t_hb_2).*max_flow+t_hb_2.*min_flow+t_hb_2.*(1-t_hb_2).*(-(max_flow-min_flow).*(1-t_hb_2)+(max_flow-min_flow).*t_hb_2);
f_hb=[f_hb_1,f_hb_2(2:end-1)];
x_hb=[x_hb_1,x_hb_2(2:end-1)];
pv = [repmat(f_hb,1,(ceil(T/cc))),f_hb_2(end)];
x_pv = 0:dt*lsf:(ceil(T/cc))*cc;
poiseuille_velocity_coefficient = interp1(x_pv,pv,(0:dt:(ceil(T/cc))*cc));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% White blood cells %%%%%%
%%%%% Space discretization %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Defining starting postition nodes on white blood cell boundary
N_wbc = 32;                     %Number of nodes on white blood cell
% Defining the shape of the white blood cells
x_initial_wbc = white_blood_cell_size*cos(((0:N_wbc-1)')*2*pi/N_wbc)';
y_initial_wbc = white_blood_cell_size*sin(((0:N_wbc-1)')*2*pi/N_wbc)'; 
% Poisson parameter for white blood cells in blood: Assigned value by other parameters.
lambda = 4/(wbc_concentration*max_flow);

% Contraction white blood cells
E_Twbc = 2*(1/lambda)*T;                    % Memory allocation: Twice the expected total white blood cells
alpha_wbc_0 = alpha * ones(1,E_Twbc);       % cell contraction forcing parameter
alpha_wbc = zeros(1,E_Twbc);
alpha_wbc(1,:) = alpha_wbc_0;
A_wbc_0 = polyarea(x_initial_wbc,y_initial_wbc);
A_wbc_min = 0.25*A_wbc_0;
A_wbc_max = 2*A_wbc_0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Virusses %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Defining area of virus population (ellipsoid)
ec = wound_center;                 %ellipsoid_center;
ea = wound_horizontal_size;        %horizontal length;
eb = wound_vertical_size;          %vertical length;
x_virus_ellipsoid = ea*cos(((0:39)')*2*pi/40)'+ec(1);
y_virus_ellipsoid = eb*sin(((0:39)')*2*pi/40)'+ec(2); 
N_virus = round(virus_density*polyarea(x_virus_ellipsoid,y_virus_ellipsoid));
%Defining population parameters
r_virus = 1;                    %Time necessary for virus reproduction 
                                %Each time interval there is a number of r_virus timesteps where reproduction is allowed 
%Define random start positions of the virusses                   
x_virus = ea*(2*rand(1,N_virus)-1)+ec(1);
y_virus = sqrt(eb^2*(1-(1/ea^2)*(x_virus-ec(1)).^2)).*(2*rand(1,N_virus)-1)+ec(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Allocating memory %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_wbc = zeros(E_Twbc,N_wbc);        %Allocation of memory for white blood cells
y_wbc = zeros(E_Twbc,N_wbc);
xc = zeros(E_Twbc,1);
yc = zeros(E_Twbc,1);
vx_wbc = zeros(E_Twbc,N_wbc); 
vy_wbc = zeros(E_Twbc,N_wbc);
wbc = 0;                            %Number of defined white blood cells
wbc_visible=(1:wbc);
distance_virus = zeros(1,N_virus);
vx_virus = zeros(1,N_virus);
vy_virus = zeros(1,N_virus);
Animation(1:K) = struct('cdata', zeros(343,436,3),'colormap',colormap(hsv(256)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Initial image and pause %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cla;  
%Plot the area of the wound
patch(x_virus_ellipsoid,y_virus_ellipsoid,(1/255)*[255, 112, 119]);
hold on
%coloring the blood vessel
patch([-L_blood_vessel/2-1,-L_blood_vessel/2-1,L_blood_vessel/2+1,L_blood_vessel/2+1],[-(thickness_blood_vessel+thickness_blood_vessel_cell)/2,(thickness_blood_vessel+thickness_blood_vessel_cell)/2,(thickness_blood_vessel+thickness_blood_vessel_cell)/2,-(thickness_blood_vessel+thickness_blood_vessel_cell)/2],1/255*[238,44,44]);
hold on;
%plotting the virusses
plot(x_virus,y_virus,'k.');   
%plotting the blood vessel cells
for i=1:N_blood_vessels
    patch(x_bvc_up(i,:),y_bvc_up(i,:),1/255*[205,16,118]);
    hold on
    patch(x_bvc_down(i,:),y_bvc_down(i,:),1/255*[205,16,118]);
    hold on
end;
%Define the final plot properties, and create an awesome animation!
axis(plot_window);               %Define x- and y-axis
whitebg([238/255,169/255,184/255]);	%Change background color of plot
%Define text in plot
text_time = text((plot_window(1)+plot_window(2))/2,4,'t = 0');
set(text_time,'FontSize',10);
set(text_time,'Color',(1/255)*[0,0,0]);
set(text_time,'FontName','Calibri');
set(text_time,'HorizontalAlignment','center');
text_N_virus = text((plot_window(1)+plot_window(2))/2,3,['Total number of virusses is ',num2str((length(x_virus))),'.']);
set(text_N_virus,'FontSize',14);
set(text_N_virus,'Color',(1/255)*[205,16,118]);
set(text_N_virus,'FontName','Calibri');
set(text_N_virus,'HorizontalAlignment','center');
text_start = text((plot_window(1)+plot_window(2))/2,2,'< Press any key to start the simulation >');
set(text_start,'FontSize',10);
set(text_start,'Color',(1/255)*[0,0,0]);
set(text_start,'FontName','Calibri');
set(text_start,'HorizontalAlignment','center');
Animation(1)=getframe;
%set(gca,'Visible','off');          %remove axes from plot
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Time integration %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:K-1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%% Image plot %%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cla;  
    %Plot the area of the wound
    patch(x_virus_ellipsoid,y_virus_ellipsoid,(1/255)*[255, 112, 119]);
    hold on
    %coloring the blood vessel
    patch([-L_blood_vessel/2-1,-L_blood_vessel/2-1,L_blood_vessel/2+1,L_blood_vessel/2+1],[-(thickness_blood_vessel+thickness_blood_vessel_cell)/2,(thickness_blood_vessel+thickness_blood_vessel_cell)/2,(thickness_blood_vessel+thickness_blood_vessel_cell)/2,-(thickness_blood_vessel+thickness_blood_vessel_cell)/2],1/255*[238,44,44]);
    hold on;
    %plotting the virusses
    plot(x_virus,y_virus,'k.');   
    %plotting the blood vessel cells
    for i=1:N_blood_vessels
        patch(x_bvc_up(i,:),y_bvc_up(i,:),1/255*[205,16,118]);
        hold on
        patch(x_bvc_down(i,:),y_bvc_down(i,:),1/255*[205,16,118]);
        hold on
    end;
    %Plot the white blood cells             
    for j = 1:wbc   
        patch(x_wbc(j,:),y_wbc(j,:),'w');   
        hold on;
        plot(xc(j,:),yc(j,:),'ko','MarkerSize',3,'MarkerFaceColor',(1/255)*[142 	56 	142]);
        hold on;
    end;
    %Define the final plot properties, and create an awesome animation!   
    axis(plot_window);                  %Define x- and y-axis
    whitebg([238/255,169/255,184/255]);	%Change background color of plot
    %Define text in plot
    text_N_virus = text((plot_window(1)+plot_window(2))/2,3,['Total number of virusses is ',num2str(length(x_virus)),'.']);
    set(text_N_virus,'FontSize',14);
    set(text_N_virus,'Color',(1/255)*[205,16,118]);
    set(text_N_virus,'FontName','Calibri');
    set(text_N_virus,'HorizontalAlignment','center');    
    text_time = text((plot_window(1)+plot_window(2))/2,4,['t = ',num2str(floor(k*dt))]);
    set(text_time,'FontSize',10);
    set(text_time,'Color',(1/255)*[0,0,0]);
    set(text_time,'FontName','Calibri');
    set(text_time,'HorizontalAlignment','center');    
    if(length(x_virus)>10*N_virus)
        text_failure = text((plot_window(1)+plot_window(2))/2,2,'The virus colony claimed victory.');
        set(text_failure,'FontSize',17);
        set(text_failure,'Color',(1/255)*[41,36,33]);
        set(text_failure,'FontName','Calibri');
        set(text_failure,'HorizontalAlignment','center');
        Animation(k+1)=getframe;
        break;
    end;        
    if isempty(x_virus)
    text_succes = text((plot_window(1)+plot_window(2))/2 ,2,'All virusses have been succesfully annihilated.');
    set(text_succes,'FontSize',17);
    set(text_succes,'Color',(1/255)*[41,36,33]);
    set(text_succes,'FontName','Calibri');
    set(text_succes,'HorizontalAlignment','center');
        Animation(k+1)=getframe;
        break;
    end;
    Animation(k+1)=getframe; 
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% Define white blood cells entering the screen  %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %New white blood cells enter following a poisson process
    if((rand<(dt/lambda)) || (k==1)) %If there comes a new white blood cell it?s location and velocity will be defined below
        wbc = wbc+1;
        % Initial position   
        xc(wbc) = (1/2)*(white_blood_cell_size-L_blood_vessel);
        yc(wbc) = (thickness_blood_vessel-white_blood_cell_size)*(betainv(rand,2,2)-(1/2)); %Using a beta(2,2)-distribution we define the spawnlocation of the white bloodcell random in the y-axis
        x_wbc(wbc,:) = xc(wbc) + x_initial_wbc;
        y_wbc(wbc,:) = yc(wbc) + y_initial_wbc;                          
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%% Define virus movement %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Annihilation by white blood cells: Annihilate virusses that are close to white blood cells
    if ~isempty(x_virus)
        distance_virus = zeros(1,length(x_virus));
        annihilation_check = zeros(1,length(x_virus));
        for j = 1:wbc 
            distance_virus(1,:) = min(((repmat(x_wbc(j,:),length(x_virus),1)-repmat(x_virus',1,N_wbc)).^2+(repmat(y_wbc(j,:),length(x_virus),1)-repmat(y_virus',1,N_wbc)).^2).^(1/2),[],2);
            annihilation_check = annihilation_check + ~isempty(x_virus).*(distance_virus(1,:) < active_range_wbc);
        end;
        x_virus(annihilation_check>0)=[];
        y_virus(annihilation_check>0)=[];
        vx_virus(annihilation_check>0)=[];
        vy_virus(annihilation_check>0)=[];    
    end;
    %Natural death: Virusses are only allowed to die if there are 2 or more.
    random_death=(length(x_virus)>1).*(rand(1,length(x_virus))<(dt/dr_virus));       
    x_virus(random_death>0)=[];
    y_virus(random_death>0)=[];
    vx_virus(random_death>0)=[];
    vy_virus(random_death>0)=[];        
    %Reproduction: Virusses can only be reproduced if there is at least 1 virus present.
    random_birth = ~isempty(x_virus).*(rem(k,1/dt)==0).*(rand(1,length(x_virus))<(dt/br_virus));
    x_virus = [x_virus,x_virus(random_birth>0)];
    y_virus = [y_virus,y_virus(random_birth>0)];
    vx_virus = [vx_virus,zeros(1,sum(random_birth>0))];
    vy_virus = [vy_virus,zeros(1,sum(random_birth>0))]; 
    %Movement of virusses is normally distributed with mean 0 and sigma^2=sigma_virus
    vx_virus = -sign(x_virus-ec(1)).*((((x_virus-ec(1))./ea).^2+((y_virus-ec(2))./eb).^2)>1)*2*sigma_virus + randn(1,length(x_virus)).*sigma_virus.^2;
    vy_virus = -sign(y_virus-ec(2)).*((((x_virus-ec(1))./ea).^2+((y_virus-ec(2))./eb).^2)>1)*2*sigma_virus + randn(1,length(x_virus)).*sigma_virus.^2;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%% Define white blood cell movement %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    wbc_visible=1:wbc;
    wbc_visible=wbc_visible(xc(1:wbc)<(0.5*L_blood_vessel+white_blood_cell_size));
    for j = wbc_visible       
        %Define new contraction parameter alpha for all blood cells
        alpha_wbc(j) = alpha_wbc_0(j)*((A_wbc_0-A_wbc_min)*(A_wbc_max-A_wbc_0))/((polyarea(x_wbc(j,:),y_wbc(j,:))-A_wbc_min)*(A_wbc_max-polyarea(x_wbc(j,:),y_wbc(j,:))));

        %Calculate virus influence on the velocity of the
        %white blood cells and add the influences using superposition.        
        vx_partial=zeros(length(x_virus),N_wbc);
        vy_partial=zeros(length(x_virus),N_wbc);        
        for s = 1:length(x_virus)
            bv_bordercheck_up = ((y_wbc(j,:)< -thickness_blood_vessel/2 & y_wbc(j,:)> -(thickness_blood_vessel_cell+(thickness_blood_vessel)/2)));
            bv_bordercheck_down = ((y_wbc(j,:)> thickness_blood_vessel/2 & y_wbc(j,:)< (thickness_blood_vessel_cell+(thickness_blood_vessel)/2)));
            bv_bordercheck = (bv_bordercheck_up|bv_bordercheck_down);
            vx_partial(s,:) = ((sum(-bv_bordercheck+1)> round(N_wbc*(4/5)))*beta+(sum(-bv_bordercheck+1)<= round(N_wbc*(4/5)))*beta/stiffness_bvc).*((1/(2*pi*D))*((x_virus(s)-x_wbc(j,:))./((x_virus(s)-x_wbc(j,:)).^2+(y_virus(s)-y_wbc(j,:)).^2)));
            vy_partial(s,:) = ((-bv_bordercheck+1)*beta+bv_bordercheck*beta/stiffness_bvc).*((1/(2*pi*D))*((y_virus(s)-y_wbc(j,:))./((x_virus(s)-x_wbc(j,:)).^2+(y_virus(s)-y_wbc(j,:)).^2))); 
        end;   
        
        distance_wbc = zeros(1,wbc);
        vx_anti_overlap_partial = zeros(1,wbc);
        vy_anti_overlap_partial = zeros(1,wbc);
        W_loop = wbc_visible;
        W_loop(W_loop==j)=[];
        for w = W_loop
            distance_wbc(w) = sqrt((xc(j)-xc(w))^2+(yc(j)-yc(w))^2);
            vx_anti_overlap_partial(w) = (1/distance_wbc(w)^0.001)*(distance_wbc(w)<1.5*white_blood_cell_size)*(1/distance_wbc(w))*(xc(j)-xc(w));
            vy_anti_overlap_partial(w) = (1/distance_wbc(w)^0.001)*(distance_wbc(w)<1.5*white_blood_cell_size)*(1/distance_wbc(w))*(yc(j)-yc(w));
        end;

        vx_wbc_virus_influence = ~isempty(x_virus)*sum(vx_partial,1);
        vx_wbc_blood_flow = -poiseuille_velocity_coefficient(k)*(y_wbc(j,:)> -(thickness_blood_vessel+thickness_blood_vessel_cell)/2 & y_wbc(j,:)< (thickness_blood_vessel+thickness_blood_vessel_cell)/2).*(y_wbc(j,:).^2-((thickness_blood_vessel+thickness_blood_vessel_cell)^2)/4);
        vx_wbc_contraction = alpha_wbc(j)*(xc(j)+x_initial_wbc-x_wbc(j,:));
        vx_wbc_random_horizontal_movement = zeros(1,N_wbc);
        vx_anti_overlap = ones(1,N_wbc)*sum(vx_anti_overlap_partial);
        vx_wbc_border = zeros(1,N_wbc);

        vy_wbc_virus_influence = ~isempty(x_virus)*sum(vy_partial,1);
        vy_wbc_blood_flow = zeros(1,N_wbc);
        vy_wbc_contraction = alpha_wbc(j)*(yc(j)+y_initial_wbc-y_wbc(j,:));
        vy_wbc_random_vertical_movement = rvm_wbc*(rand-0.5)*(sum((y_wbc(j,:)> -thickness_blood_vessel/2 & y_wbc(j,:)< thickness_blood_vessel/2))==N_wbc);
        vy_anti_overlap = ones(1,N_wbc)*sum(vy_anti_overlap_partial);
        vy_wbc_border = (1/500)*sum(-bv_bordercheck+1) * stiffness_bvc * (-max(bv_bordercheck_down>0)+max(bv_bordercheck_up>0));
        
        vx_wbc(j,:) = vx_wbc_virus_influence + vx_wbc_blood_flow + vx_wbc_contraction + vx_wbc_random_horizontal_movement + vx_anti_overlap + vx_wbc_border;
        vy_wbc(j,:) = vy_wbc_virus_influence + vy_wbc_blood_flow + vy_wbc_contraction + vy_wbc_random_vertical_movement + vy_anti_overlap + vy_wbc_border;
    end; 
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%% Euler forward %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    x_virus = x_virus + ~isempty(length(x_virus)).*dt*vx_virus;
    y_virus = y_virus + ~isempty(length(x_virus)).*dt*vy_virus;    
    x_wbc = x_wbc + dt*vx_wbc;
    y_wbc = y_wbc + dt*vy_wbc;    
    xc = (1/N_wbc)*sum(x_wbc,2);      %Define cell center
    yc = (1/N_wbc)*sum(y_wbc,2);
end;

Animation(k+2:K) = [];
%Save movie as .avi file.
%movie2avi(Animation, 'filename.avi','fps',80) 
%Compress file with virtualdub

 