Team:Exeter/Modelling BM01
<!DOCTYPE html
PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
Contents
- BM_01_09
- Setting our deafault parameters and variables
- The coordinate generation and the main loop
- Generation of the inital starting points, calls startpoint.
- Coordinate generation and confinement checking
- Joining and splitting
- Bouncing
- Startpoint
- Checking function
- Joining function
- Splitter function
- Plotting function
- Bouncing function
- Counting function
- Code needed to produce a .gif file of the simulation output
- End Of The Code
BM_01_09
Our final code for the simulation of our cell free kit. This is after the restructure implimented after meetings with both Jonathan Fieldsend and the lab team. There are still improvements and optimisations to be made these are discussed below. More information can be found at https://2015.igem.org/Team:Exeter/Modeling. The function take inputs of: * t - Number of toeholds * r - Number of RNA's * N - Number of time steps * T - temperature -> this is the parameter scanning variable It outputs GFPcount, this depends on the parameter scanning variable choosen.
function [GFPcount] = BM_01_09(t,r,N,T)
Setting our deafault parameters and variables
These are the parameters used for the basic setup of Brownian motion as well as the contianment to a tube. Containment changed to a cylinder to emulate the lab, they are using a well plate now. All units are SI units unless otherwise stated.
rng('shuffle'); eta = 1.0e-3; % viscosity of water in SI units (Pascal-seconds) kB = 1.38e-23; % Boltzmann constant %T = 293; % Temperature in degrees Kelvin tau = .1; % time interval in seconds d_t=5.1e-8; % diameter in meters of toehold d_r=1.7e-8; % of RNA d_c=5.1e-8; % of toehold-RNA complex D_t = kB * T / (3 * pi * eta * d_t); %diffusion coefficient toehold D_r = kB * T / (3 * pi * eta * d_r); %diffusion coefficient RNA D_c = kB * T / (3 * pi * eta * d_c); %diffusion coefficient complex p_t = sqrt(2*D_t*tau); p_r = sqrt(2*D_r*tau); p_c = sqrt(2*D_c*tau); %CONFINEMENT - height can be changed depending on the volume of %the solution (rather than the total possible volume) A = 2.5e-2; %binding distance, real world value 2.5e-10 cylinder_radius=3.34e-3; tube_height_max=5.45e-3; tube_height_min=-5.45e-3; cone_height=18e-3; %GIF stuff theta = 0; % changes the viewing angle in the GIF %Figure Stuff -> setting up the grid for the simulation % figure() % %axis([-0.00005 0.00005 -0.00005 0.00005 -0.00005 0.00005]); % axis([-5e-3 5e-3 -5e-3 5e-3 -8e-3 8e-3]) %changed % grid on % grid MINOR % set(gcf, 'Position', [100 10 600 600]) % xlabel('X-axis') % ylabel('Y-axis') % zlabel('Z-axis') % hold on %Choosing the number of complexs based on the lower of t and r if t>=r c=r; else c=t; end %Checking variables and initalising to zeros blank=[0,0,0]; points={}; joinstatus=zeros(N,t+r); %Initalising points to zeros for j=1:N for k=1:t+r+(2*c) points{j,k}=blank; end for k=t+r+2:2:t+r+(2*c) points{j,k}=0; end end bouncepoints=points;
The coordinate generation and the main loop
The coordinates are now generated on the fly at every time step, this prevents the need to pass around large matrices. This is also the main loop at which every required function is called at each time step.
%The main for loop -> loops over all the time steps for j=1:N
Generation of the inital starting points, calls startpoint.
if j==1 [coords, startposition] = startpoint(); check=zeros(c,3); else
Coordinate generation and confinement checking
The coordinates are checked against the height of the tube frist then checkxy is called to check whether it has left the cyclinder.
%Toehold for i=1:t if any(any(check==i))~=1 for k=1:3 coords(k+3,i)=coords(k,i); coords(k,i)=(p_t*randn(1))+coords(k+3,i); end % checking against the height of the tube. if coords(3,i)>=tube_height_max coords(3,i)=tube_height_max; elseif coords(3,i)<=tube_height_min coords(3,i)=tube_height_min; end %calling checkxy to check x and y [coords,bouncepoints]=checkxy(cylinder_radius, ... coords, i, bouncepoints); end end %RNA for i=t+1:t+r if any(any(check==i))~=1 for k=1:3 coords(k+3,i)=coords(k,i); coords(k,i)=(p_r*randn(1))+coords(k+3,i); end % checking against the height of the tube. if coords(3,i)>=tube_height_max coords(3,i)=tube_height_max; end if coords(3,i)<=tube_height_min coords(3,i)=tube_height_min; end %calling checkxy to check x and y [coords,bouncepoints]=checkxy(cylinder_radius, ... coords, i, bouncepoints); end end %Complex for i=t+r+1:t+r+c if check(i-(t+r),1)~=0 && check(i-(t+r),2)~=0 for k=1:3 coords(k+3,i)=coords(k,i); coords(k,i)=(p_c*randn(1))+coords(k+3,i); end % checking against the height of the tube. if coords(3,i)>=tube_height_max coords(3,i)=tube_height_max; end if coords(3,i)<=tube_height_min coords(3,i)=tube_height_min; end %calling checkxy to check x and y [coords,bouncepoints]=checkxy(cylinder_radius, ... coords, i, bouncepoints); end end
end
Joining and splitting
This is were the joining and splitting functions are called, check is used to determine whether a complex is formed or not and hence which function is called.
for q=1:c if check(q,3)==0 && check(q,1)==0 && check(q,2)==0 && j~=1 %variable to make sure joiner and splitter dont happen in same time step jointime=j; [coords,check, startposition, points] = joiner(coords, ... check, startposition, points); elseif check(q,3)~=0 check(q,3)=check(q,3)+1; end if check(q,1)~=0 && check(q,2)~=0 && j~=jointime [coords, check, startposition, points] = splitter(q, ... coords, check, startposition, points); end end
Bouncing
for f=1:t+r+c if f<t+r+1 if bouncepoints{j,f}==0 %not a bouncer matcoords=[coords(1,f),coords(2,f),coords(3,f)]; points{j,f}=matcoords; %couldn't have these statements on the same line for some reason elseif bouncepoints{j,f}~=0 if size(bouncepoints{j,f},2)==3 matcoords=[coords(1,f),coords(2,f),coords(3,f) ... bouncepoints{j,f}(1,1), ... bouncepoints{j,f}(1,2),bouncepoints{j,f}(1,3)]; points{j,f}=matcoords; end end end if f>=(t+r+1) if f==(t+r+1) g=f; h=f+1; end if bouncepoints{j,g}==0 %not a bouncer matcoords=[coords(1,f),coords(2,f),coords(3,f)]; points{j,g}=matcoords; elseif bouncepoints{j,g}~=0 if size(bouncepoints{j,f},2)==3 matcoords=[coords(1,f),coords(2,f),coords(3,f) ... bouncepoints{j,g}(1,1), ... bouncepoints{j,g}(1,2),bouncepoints{j,g}(1,3)]; points{j,g}=matcoords; end end g=g+2; if j>1 if points{j-1,h}(1,1)~=0 if points{j,h}==0 points{j,h}=1; end end end h=g+1; end end if j==N counter(points); finish=1; end
end
Startpoint
Creates a start point definitely inside the dimensions of container, this is randomly generated. A vector called startpositon is made so the particles are located around this.
function [coords, startposition] = startpoint() %each column is a toehold, with six rows, %for current xyz and previous xyz coords=zeros(6,(t+r+c)); for m=1:t+r coords(3,m)=(tube_height_min)+ ... ((tube_height_max-tube_height_min)*rand(1)); if coords(3,m)>=0 coords(1,m)=(-cylinder_radius)+ ... ((cylinder_radius-(-cylinder_radius))*rand(1)); coords(2,m)=(-cylinder_radius)+ ... ((cylinder_radius-(-cylinder_radius))*rand(1)); end if coords(3,m)<0 cone_radius=(0.26*coords(3,m))+4.6e-3; coords(1,m)=(-cone_radius)+ ... ((cone_radius-(-cone_radius))*rand(1)); coords(2,m)=(-cone_radius)+ ... ((cone_radius-(-cone_radius))*rand(1)); end end startposition=coords; end
Checking function
Checks the coordinates are within the boundaries of the eppendorf
%Function that checks whether the particle is inside of the tube %for its calculated z-coordinate at the point of contact in the %tube (in cone: via line equation intersects of path and boundary %// in cylinder: line equation of tube boundary). %From this point of contact, the new X and Y coordinates are %calculated for the "exit point" and then subsequently, the new %resultant xyz can be calculated function [coords,bouncepoints]=checkxy(radius,coords,i,bouncepoints) if (coords(1,i)^2)+(coords(2,i)^2)>=(radius^2) %red box in (Maths for exitx.png explains derivation) %setting exitX/Y at boundary of cylinder grad=abs((sqrt((coords(1,i)^2)+(coords(2,i)^2))- ... sqrt((coords(4,i)^2)+(coords(5,i)^2)))/(coords(3,i) ... -coords(6,i))); %confirmed to be correct if coords(1,i)<0 exitX=-sqrt((radius^2)/((grad^2)+1)); else exitX=sqrt((radius^2)/((grad^2)+1)); end if coords(2,i)<0 exitY=-(grad*sqrt(((radius^2)/((grad^2)+1)))); else exitY=grad*sqrt(((radius^2)/((grad^2)+1))); end Px=coords(1,i); Py=coords(2,i); lastx=coords(4,i); lasty=coords(5,i); trajectory_gradient=(Py-lasty)/(Px-lastx); %m1 tangent_gradient=(-exitX/exitY); %m2 theta_bounce=atan(trajectory_gradient); phi_2=atan(tangent_gradient); phi_1=theta_bounce-phi_2; exitPDist=sqrt(((Px-lastx)^2)+((Py-lasty)^2)); G_length=exitPDist*cos(phi_1); Gx=exitX+(G_length*cos(phi_1)); Gy=exitY+(G_length*cos(phi_1)); Cx=Px-Gx; Cy=Py-Gy; newX=Px-(2*Cx); newY=Py-(2*Cy); perpendicular_gradient=(exitY/exitX); %z-intercept with boundary: (see book) prad=sqrt((coords(1,i)^2)+(coords(2,i)^2)); lastrad=sqrt((coords(4,i)^2)+(coords(5,i)^2)); Dr=prad-lastrad; pZ=coords(3,i); lastZ=coords(6,i); Dz=pZ-lastZ; Gr=radius-lastrad; Gz=Gr*tan(acos(Dr/sqrt((Dr^2)+(Dz^2)))); if pZ<lastZ exitZ=lastZ-Gz; else exitZ=lastZ+Gz; end %write new points directly into points in a way that joiner and %splitter can still work if i<=t+r bouncepoints{j,i}=[exitX exitY exitZ]; elseif i>t+r bouncepoints{j,(2*i)-1}=[exitX exitY exitZ]; end %update for next point coords(1,i)=newX; coords(2,i)=newY; coords(3,i)=coords(3,i); end %What if the particle attempts to fall out of the bottom of the tube? end
Joining function
This function has a threshold for joining, a joining probability. If this is met the toehold and RNA's are joined. The check vector changes to indicate this, the toeholds and RNA's are spotted being made and a complex line is made instead. startingposition is also updated to indicate the new starting point of the complex line.
function [coords, check, startposition, points] = joiner(coords, ... check, startposition, points) colshift=0; joinprob = randi([1 10],1); if joinprob >=3 for n=1:t for m=t+1:r+t if any(any(check==n))==1 || any(any(check==m))==1 continue else % if ((((tx(j,k)-rx(j,m))^2)+((ty(j,k)-ry(j,m))^2)+((tz(j,k)-rz(j,m))^2))<=(A^2) || (check_r(1,m)~=0 && check_t(1,k)~=0)) && (j~=1) && delay(1,n)==0 %dist2 on sub matrices to find distances between for joining if (((coords(1,n)-coords(1,m))^2)+((coords(2,n)-... coords(2,m))^2)+((coords(3,n)- ... coords(3,m))^2))<=(A^2) for p=1:c if check(p,1)~=0 && check(p,2)~=0 continue else if check(p,1)==0 && check(p,2)==0 coords(1,t+r+p)=(coords(1,n)+ ... coords(1,m))/2; coords(2,t+r+p)=(coords(2,n)+ ... coords(2,m))/2; coords(3,t+r+p)=(coords(3,n)+ ... coords(3,m))/2; startposition(1,t+r+p)= ... coords(1,t+r+p); startposition(2,t+r+p)=... coords(2,t+r+p); startposition(3,t+r+p)=... coords(3,t+r+p); check(p,1)=n; check(p,2)=m; if colshift==0; colshift=1; end points{j,t+r+p+colshift}=[1,n,m]; colshift=colshift+1; joinstatus(j,n)=1; joinstatus(j,m)=1; end break end end end end end end end end
Splitter function
Splitter is similar to joiner but the opposite happnens. If a threshold is reached the complex will split back into a toehold and RNA. These lines are now produced again, check is updated to reflect this and finally startposition is changed.
function [coords, check, startposition, points] = ... splitter(q, coords, check, startposition, points) split = randi([1 10],1); if split<=4 if check(q,1)~=0 && check(q,2)~=0 toehold=check(q,1); rna=check(q,2); startposition(1,toehold)=coords(1,t+r+q); startposition(2,toehold)=coords(2,t+r+q); startposition(3,toehold)=coords(3,t+r+q); startposition(1,rna)=coords(1,t+r+q); startposition(2,rna)=coords(2,t+r+q); startposition(3,rna)=coords(3,t+r+q); coords(1,rna)=coords(1,t+r+q); coords(2,rna)=coords(2,t+r+q); coords(3,rna)=coords(3,t+r+q); coords(1,toehold)=coords(1,t+r+q); coords(2,toehold)=coords(2,t+r+q); coords(3,toehold)=coords(3,t+r+q); points{j,t+r+(2*q)}=[0,toehold,rna]; check(q,1)=0; check(q,2)=0; check(q,3)=-5; end end end
Plotting function
function [] = plotter(points) %Conditions from plotting drawn from points array for row=1:N if row==1 %startpoints for col=1:t+r plot3(points{row,col}(1,1), points{row,col}(1,2), ... points{row,col}(1,3),'kx') end else for altcol=t+r+2:2:t+r+(2*c) if points{row,altcol}(1,1)==1 && size(points{ ... row,altcol},2)==3 && row>1 if points{row-1,altcol}(1,1)==0 %definitely joinpoint %plot blue using row-1 to row to joinpoint column, column selected from points{row,col}(1,2)) plot3([points{row-1,(points{row,altcol}(1,2))}(1,1),... points{row,altcol-1}(1,1)], ... [points{row-1,(points{row,altcol}(1,2))}(1,2),... points{row,altcol-1}(1,2)],... [points{row-1,(points{row,altcol}(1,2))}(1,3),... points{row,altcol-1}(1,3)], 'b') %plot red using row-1 to row to joinpoint column, column selected from points{row,col}(1,2) plot3([points{row-1,(points{row,altcol}(1,3))}(1,1),... points{row,altcol-1}(1,1)], ... [points{row-1,(points{row,altcol}(1,3))}(1,2),... points{row,altcol-1}(1,2)],... [points{row-1,(points{row,altcol}(1,3))}(1,3), ... points{row,altcol-1}(1,3)], 'r') %plot circle at joinpoint plot3(points{row,altcol-1}(1,1), ... points{row,altcol-1}(1,2), ... points{row,altcol-1}(1,3), 'ko') end if points{row-1,altcol}(1,1)==1 %defintely splitpoint %do nothing end elseif (points{row,altcol}(1,1)==1 && ... size(points{row,altcol},2)==1) || ... (points{row,altcol}(1,1)==0 && ... size(points{row,altcol},2)==3) %plot green from row-1 to row if size(points{row-1,altcol-1},2)==3 plot3([points{row-1,altcol-1}(1,1),... points{row, altcol-1}(1,1)], ... [points{row-1,altcol-1}(1,2), ... points{row, altcol-1}(1,2)], ... [points{row-1,altcol-1}(1,3),... points{row, altcol-1}(1,3)],'g') elseif size(points{row,altcol-1},2)==6 %Bouncing of the side points index for this is 6 numbers, the xyzpoint before hitting and the xyzpoint on the boundary bounceplot(points, row, altcol-1, 'g') end end if points{row,altcol}(1,1)==0 %currently split if size(points{row-1,altcol},2)==3 &&... points{row-1,altcol}(1,1)==0 %prevents incorrect plotting in the case of join/split in consecutive timesteps %plot from split column to blue column plot3([points{row-1,altcol-1}(1,1),... points{row,(points{row-1,altcol}(1,2))}(1,1)],... [points{row-1,altcol-1}(1,2), ... points{row,(points{row-1,altcol}(1,2))}(1,2)],... [points{row-1,altcol-1}(1,3), ... points{row,(points{row-1,altcol}(1,2))}(1,3)],... 'b') %plot from split column to red column plot3([points{row-1,altcol-1}(1,1), ... points{row,(points{row-1,altcol}(1,3))}(1,1)],... [points{row-1,altcol-1}(1,2), ... points{row,(points{row-1,altcol}(1,3))}(1,2)],... [points{row-1,altcol-1}(1,3),... points{row,(points{row-1,altcol}(1,3))}(1,3)],... 'r') %plot star at splitpoint plot3(points{row-1,altcol-1}(1,1), points... {row-1,altcol-1}(1,2), points... {row-1,altcol-1}(1,3), 'k*') %plot continue %might not be necessary end end for col=1:t+r if points{row,col}(1,1)~=points{row-1,col}(1,1) ... && points{row,col}(1,2)~=points{row-1,col}... (1,2) && points{row,col}(1,3)~=points{row-1,... col}(1,3) % have to be specific due to 1x6 object when bouncing occurs if size(points{row,altcol},2)~=3 %determine column in t or r for colour of line if col<=t %plot blue for toeholds if ((size(points{row-1,altcol},2)==3 ... && points{row-1,altcol}(1,1)==0)... && points{row-1,altcol}(1,2)==col) %if just split, don't plot from row-1 to row in t column since plotting from split point has just occurred continue elseif points{row,altcol}==0 && size... (points{row-1,altcol},2)~=3 if size(points{row-1,col},2)==3 plot3([points{row-1,col}(1,1),... points{row,col}(1,1)],... [points{row-1,col}(1,2),... points{row,col}(1,2)],... [points{row-1,col}(1,3),... points{row,col}(1,3)],'b') elseif size(points{row,col},2)==6 %Bouncing of the side points index for this is 6 numbers, the xyzpoint before hitting and the xyzpoint on the boundary bounceplot(points, row, col, 'b') end end end if col>t %plot red for rna(trigger) if ((size(points{row-1,altcol},2)==3 ... && points{row-1,altcol}(1,1)==0) ... && points{row-1,altcol}(1,3)==col) %if just split, don't plot from row-1 to row in r column since plotting from split point has just occurred continue elseif points{row,altcol}==0 && size ... (points{row-1,altcol},2)~=3 if size(points{row-1,col},2)==3 plot3([points{row-1,col}(1,1),... points{row,col}(1,1)],... [points{row-1,col}(1,2),... points{row,col}(1,2)],... [points{row-1,col}(1,3),... points{row,col}(1,3)],'r') elseif size(points{row,col},2)==6 %Bouncing of the side points index for this is 6 numbers, the xyzpoint before hitting and the xyzpoint on the boundary bounceplot(points, row, col, 'r') end end end end end end end end end end
Bouncing function
function [points] = bounceplot(points, row, col, style) plot3([points{row-1,col}(1,1) points{row,col}(1,4)], ... [points{row-1,col}(1,2) points{row,col}(1,5)], ... [points{row-1,col}(1,3) points{row,col}(1,6)],style) plot3([points{row,col}(1,4) points{row,col}(1,1)], ... [points{row,col}(1,5) points{row,col}(1,2)], ... [points{row,col}(1,6) points{row,col}(1,3)],style) end
Counting function
This function is related to the generation of a GFP output of our system, it is used in combination with the parameter scanning aspect of our simulation.
function counter(points) %simplfied status of each toehold and trigger drawn from points array for jstatcol=1:t+r for jstatrow=2:N if points{jstatrow,jstatcol}(1)==points{jstatrow-1, ... jstatcol}(1) && points{jstatrow,jstatcol}(2)== ... points{jstatrow-1,jstatcol}(2) && points{ ... jstatrow,jstatcol}(3)==points{jstatrow-1,jstatcol}(3) joinstatus(jstatrow,jstatcol)=1; % if joinstatus{jstatrow-1,jstatcol}==0 % joinstatus{jstatrow-1,jstatcol}=1; % end else if joinstatus(jstatrow,jstatcol)~=1 joinstatus(jstatrow,jstatcol)=0; end if points{jstatrow-1,jstatcol}==1 joinstatus(jstatrow,jstatcol)='split'; end end end end totalgreen=0; for col=1:t totalgreen=totalgreen+sum(joinstatus(:,col)); end totalgreentime=totalgreen*tau; GFPrate=20; GFPcount=totalgreentime*GFPrate; %GFPconc=GFPcount/eppendorfvolume; end
Code needed to produce a .gif file of the simulation output
function jiff(row) change = 360/N; % the size of the angle change % gif utilities set(gcf,'color','w'); % set figure background to white drawnow; frame = getframe(gcf); im = frame2im(frame); [imind,cm] = rgb2ind(im,256); outfile = '17_07_v1.gif'; % adjusting the viewing the angle view(theta,45); theta = theta + change; % On the first loop, create the file. In subsequent loops, append. if row==2 imwrite(imind,cm,outfile,'gif','DelayTime',0,'loopcount',inf); else imwrite(imind,cm,outfile,'gif','DelayTime',0,'writemode', ... 'append'); end end
end
End Of The Code
This is the final simulation developed by the Univeristy Of Exeter iGEM team 2015. Developed mainly by Amy, Dan and Todd.