% 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);