% MATLAB code for the article "Macromolecular crowding directs the motion
% of small molecules inside cells" by Stephen Smith and Ramon Grima.
%
% This code simulates the Brownian motion of a single sphere of radius r in
% an arbitrary distribution phi of spheres of radius R. The initial
% distribution phi is specified by the user, and random non-overlapping
% sphere locations are picked to fit the distribution. The Cichocki-Hinsen
% algorithm is used for the particle motion. The output is a vector Ystore
% which contains the trajectory of the small particle.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameter specification
R=0.05; %The radius of the large particles
epsilon=0.1; %The ratio of particle sizes r/R
r=epsilon*R; %The radius of the small particle
D=1; %Small particle diffusion coefficient
DC=D*epsilon; %Large particle diffusion coefficient
dt=1e-05; %Simulation time-step, should be as small as feasible
nSim=1e02; %Number of simulation time steps.
x=-10:(2*R):10; %The grid used to define the inital particle positions
phi=0.1*normpdf(x,0,1)/max(normpdf(x,0,1)); %The initial large particle
% density. phi represents the local proportion of volume occupied, as a
% function of position.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This section uses the grid x to place the large particles according to
% the distribution phi
l=1;
X=[]; %Array containing particle coordinates
for i=1:length(x) %Iterate over the grid points
n=round(2*R*phi(i)/((4/3)*pi*R^3)); %Number of large particles in gridpoint i
for j=1:n %Iterate over the particles at the gridpoint
prop=[x(i)+(0.1*rand-0.05);R+(1-2*R)*rand;R+(1-2*R)*rand]; %proposed large particle position.
if l>1
[~,dx]=knnsearch(X',prop');
while dx<(2*R) %Keep choosing new positions until there are no intersections with other particles
prop=[x(i)+(0.1*rand-0.05);R+(1-2*R)*rand;R+(1-2*R)*rand];
[~,dx]=knnsearch(X',prop');
end
end
X=[X,prop];
l=2;
end
end
NC=length(X); %Total number of large particles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This section places the small particle at 0 on the x axis, but not
% overlapping any large particles.
Y=[0;r+(1-2*r)*rand;r+(1-2*r)*rand];
[~,dx]=knnsearch(X',Y');
while dx<(r+R)
Y=[0;r+(1-2*r)*rand;r+(1-2*r)*rand];
[~,dx]=knnsearch(X',Y');
end
YStore=Y(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This section simulates the diffusion of the particles. At each time step,
% a new position is proposed for each particle in turn. If the proposed
% position would result in a collision, the move is rejected, otherwise it
% is accepted.
h=waitbar(0);
for i=1:nSim %Iterate over each time step
waitbar(i/nSim,h,[num2str(floor(100*i/nSim)),'% complete...']);
prop=normrnd(0,sqrt(2*D*dt),[3,1]); %Proposed small particle position
[~,dd]=knnsearch(X',(Y+prop)');
if Y(2)+prop(2)<(1-r)&&Y(2)+prop(2)>r&&Y(3)+prop(3)<(1-r)&&Y(3)+prop(3)>r&&dd>(r+R)
Y=Y+prop; %Update small particle position
end
for j=1:NC %Iterate over each large particle
prop=normrnd(0,sqrt(2*DC*dt),[3,1]); %Proposed large particle position
[~,dx]=knnsearch(X(:,[1:(j-1),(j+1):NC])',(X(:,j)+prop)');
[~,dy]=knnsearch(Y',(X(:,j)+prop)');
if X(2,j)+prop(2)<(1-R)&&X(2,j)+prop(2)>R&&X(3,j)+prop(3)<(1-R)&&X(3,j)+prop(3)>R&&dx>(2*R)&&dy>(r+R)
X(:,j)=X(:,j)+prop; %Update large particle position
end
end
Ystore(i+1)=Y(1); %Store the x-coordinate of the small particle.
end
close(h);