   % this is a simple one dimensional agent-based simulation of diffusion via the run and
   % tumble process
% the code can be modified to add multiple states and state transitions.
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);
clear
N = 2000; % number of particles
Nstep = 2000;  %This is the total number of time steps to take
dt = 0.025;  % time step size 
v =[-1;0;1]; % the ballistic velocity of the particles in different states 
kon = 1; % assign change rate values
koff = 1;
alpha = 0.5; % the probability of going left rather than right
t=0;
x= zeros(N,1);  %the particle position
s = randi(3,[N,1]); %the initial state of each particle
% the three states are leftward motion s=1, resting (no motion) s=2,
% rightward motion s=3.
pltflag = 0;  % make a plot of the time course or not
  % Beware:
  %pltflag =1 makes the simulation MUCH slower!
% the change rate vector
crate = [koff,kon,koff];  % these are the rates for changing OUT of each state
c = [0;1;0]; %these are used to determine the new states after a change
for j = 1:Nstep
    % first move the particles
    x = x +v(s)*dt;
    t = t+dt;
% now decide if states should change
 % pick N random numbers
 R = rand(N,1);
del = find(crate(s)'*dt>=R);  % find the particles that should change state
ndel = length(del);
s(del) = c(s(del)).*sign(rand(ndel,1)-alpha) +2; % assign the new state value
%this completes one cycle

%make a plot of current positions
if(pltflag==1)
    x2 = find(s==2);
x1 = find(s==1);
x3 = find(s==3);
figure(1)
plot(t*ones(length(x1)),x(x1),'b*')
hold on
plot(t*ones(length(x2)),x(x2),'r*')
plot(t*ones(length(x3)),x(x3),'k*')
xlabel('time','fontsize',20)
ylabel('x','fontsize',20)
end
% keep track of the mean squared displacement
T(j) = t;
Xsq(j) = sum(x.^2)/N;
end

hold off

% make a histogram of the positions
figure(2)
Mx = max(abs(x));
dx = 1/20;
binc = [-1:dx:1]*Mx;
    hist(x,binc)
    xlabel('x','fontsize',20)
    % plot the mean squared displacement
    figure(3) 
    plot(T,Xsq,'*')
    xlabel('Time','fontsize',20)
    ylabel('Mean Squared Displacement','fontsize',20)
