Difference between revisions of "Team:Exeter/Modelling BM01"

Contents

BM_07_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_07_09_insilicobinding(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
    D_r = kB * T / (3 * pi * eta * d_r);
    D_c = kB * T / (3 * pi * eta * d_c);
    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 of the eppendorf)
    A = (3.5e-10)*2; %binding distance, default at 1e-7, real world value .5e-10

    cylinder_radius=3.34e-6; %radius of F-Well plate well in metres (3.34e-3metres) %changed
    tube_height_max=7.13e-7; %height of F-Well plate well to fill i billionth of 50microlitres (1.426e-3metres)
    tube_height_min=-7.13e-7;
    %Changeable to suit the container of your system
    %Total height is split in half, with one half above the positive xy plane and one half below
    cone_height=18e-3; %unused


    %GIF stuff
    theta = 0; % changes the viewing angle

    %Figure Stuff
%     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();
%             c_joined=zeros(1,c);
%             c_split=zeros(1,c);
            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
        for q=1:c
            if check(q,3)==0 && check(q,1)==0 && check(q,2)==0 && j~=1
                jointime=j; %variable to make sure joiner and splitter dont happen in same time step
                [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 %couldn't have these statements on the same line for some reason
                        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
            plotter(points);
            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()
        coords=zeros(6,(t+r+c)); %each column is a toehold, with six rows, for current xyz and previous xyz
%         startposition=zeros(6,(t+r+c));
%         coords=mat2dataset(coords);
        for m=1:t+r
            coords(3,m)=(tube_height_min)+((tube_height_max-tube_height_min)*rand(1));
            coords(1,m)=(-cylinder_radius)+((cylinder_radius-(-cylinder_radius))*rand(1));
            coords(2,m)=(-cylinder_radius)+((cylinder_radius-(-cylinder_radius))*rand(1));

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


            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);
%                 shifted_perpendicular_line=(perpendicular_gradient*EX)-((exitY+Px)/exitX)+Py;

                %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
    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;
        %Joining probability calculated from in silico data of free energy of complex from NUPACK and
        %equation of polynomial fit to normalised curve gives probability
        %of binding (or rather gives the threshold to enable a successful
        %binding event which a randomly generated number is tested against)
        joinevent = (randi([1 10000000],1))/10000000;
        Tempcorrection = T-273;
        jointhreshold = ((4e-07)*(Tempcorrection^4)) - ((6e-05)*(Tempcorrection^3)) + (0.0023*(Tempcorrection^2)) - (0.0072*Tempcorrection) + 0.3745;
        if joinevent >= jointhreshold
            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
                        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