Instructions to run simulation for crawling nematode sperm cell:

To run the program, you need to download all these files into a Matlab folder in your "Current Directory." The main program is “carcasse.m”, which calls all other subroutines as functions.”

Then type in "carcasse" in the prompt.

It will then ask you how much time for the simulation. Enter a value between .4 and 2 (the interval (.4,2) ). You can probably choose a higher time value than 2, but it might take longer to compile

It will then ask you how many filaments you want to have. Enter an odd integer greater than 7. Note it must be ODD.

It will then ask you to choose between an ellipse with equal angles or equal distances. We are uncertain if both work, however, I have always used value 1, using equal angles (it is what we used for the presentation).

It will then ask you to enter the coefficient which will yield the shape of the ellipse (initial condition). Enter a value greater than 0 and less than 1 (the interval (0,1) ).

It will then ask you how many times you would like to turn. I believe you can only enter a 0 (for no turns) or a 1 (for 1 turn).

It will then ask for the angle of the turn (dictated by the external pH). I have used values like "pi/30" or "pi/60." You can play around with the angle size, but it must be in radian.

The program should then run and you should be able to watch the video.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%The carcasse of the program...SAVE AS “carcasse.m”

% carcasse.m

clear;

%variables THESE ARE COORDINATES OF TOP AND BOTTOM OF CELL%%%

tau = 0.01;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%

%%%%%%%%%%%%%%%%INPUTING DATA%%%%%%%%%%%%%%%%%%%%%%%

time = input('\nEnter the ammount of time for simulation\n');

fil = input('\nEnter the number of filaments\n');

fprintf('\nChoices you have: \n');

fprintf('\nElipse, equal angles: Press 1\n');

fprintf('\nElipse, equal y distances: Press 2\n');

choice = input('\nYour choice:\n');

if (choice~=1 && choice~=2 )

error = ('Error: Choice ~= 1/2');

end

%%%the rightmost point of the initial elipse

coeff = input('\nChoose the coefficient\n');

lam_size = coeff*1.1;

x0(1:time/tau) = 0;

y0(1:time/tau) = 1;

xf(1:time/tau) = 0;

yf(1:time/tau) = -1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%INITIALIZE VECTOR OF TURNING ANGLES%%%%%%%%%%%%

turn_count = input('\nHow many times would you like to turn?\n');

if (turn_count>0)

turn_angles(1:turn_count) = zeros;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%ASKS FOR THE INPUT OF DEGREES OF TURN%%%%%%%%%%

if (turn_count>0)

for i =1:turn_count

turn_angles(i) = input('\nInput an angle:\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%CHECKING THAT SUM, ANGLE IS B/W +-pi/2%%%%%%%

for i = 1:turn_count

if (sum(turn_angles(1:turn_count))>pi/2 || sum(turn_angles(1:turn_count))<-pi/2 ...

|| turn_angles(i)>pi/2 || turn_angles(i)<-pi/2)

error = ('\nError w/ the angles\n');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%INITIALIZE THE X, Y COORDINATE MATRIX%%%%%%

nnx(fil, time/tau) = zeros;

nny(fil, time/tau) = zeros;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%INITIALIZE THE SHAPE%%%%%%%%%%%%%%%%%

[nnx(:, 1), nny(:, 1)]= shape_init(choice, coeff, fil);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%CALCULATE THE TIME CB IS TRIGERED%%%%%%%%%%

start = start_move(lam_size, tau);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%INITIALIZING THE TURNING TIMES%%%%%%%%%%%%%%%%%

turn_times(1:turn_count) = zeros;

if (turn_count>0)

for i = 1:turn_count

turn_times(i) = input('Input a timestep');

if (i>1 && turn_times(i)<turn_times(i-1))

error = ('Error, Time steps have to be increasing');

end

if (turn_count>0)

if (turn_times(1)<start)

error = ('Error, First time step needs to be after the cell body is trigered');

end

if (turn_times(turn_count)>time/tau)

error = ('Error,Last step needs to be before the body stops moving');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('# of start step = %f', start/tau);

fprintf('# of time steps = %f', time/tau);

%%%%%%%%%%%MOVEMENT W/O CB%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% here we assume that the ends of the cell stay in the same place

for i = 2:min((start/tau-1), time/tau)

[add, del] = check_for_add(x0(i-1), y0(i-1), xf(i-1), yf(i-1)...

, nnx(:, i-1), nny(:, i-1), fil, 0); %add can equal to 0 through fil-1

%0-nothing is added, add>0- filament is added between add and add +1

if add == -1

[nnx(:, i), nny(:, i)]= move_wob(nnx(:, i-1), nny(:, i-1), fil, tau, lam_size);

else

[nnx(:, i), nny(:, i)]= add_fil(nnx(:, i-1), nny(:, i-1)...

, x0(i), y0(i), xf(i), yf(i), add, del);

fprintf('\nI have added a filament number %f after %f timestep\n',add, i);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%REST OF MOVEMENT IF NO TURNS%%%%%%%%%%%%%%%

if turn_count ==0

for i = start/tau:time/tau

[add, del] = check_for_add(x0(i-1), y0(i-1), xf(i-1), yf(i-1),...

nnx(:, i-1), nny(:, i-1), fil, 0); %add can equal to 0 through fil-1

%0-nothing is added, add>0- filament is added between add and add +1

if add == -1

[nnx(:, i), nny(:, i)]= move_wb(0, nnx(:, i-1),...

nny(:, i-1), fil, tau, lam_size, x0(i-1), y0(i-1), xf(i-1), yf(i-1));

edx = sum(nnx(:, i)-nnx(:, i-1))/fil;

x0(i) = x0(i-1)+tau*.9;

y0(i) = 1;

xf(i) = xf(i-1)+tau*.9;

yf(i) = -1;

else

[nnx(:, i), nny(:, i)]= add_fil(nnx(:, i-1), nny(:, i-1),...

x0(i-1), y0(i-1), xf(i-1), yf(i-1), add, del);

fprintf('\nI have added a filament number %f after %f timestep\n',add, i);

x0(i) = x0(i-1);

xf(i) = xf(i-1);

y0(i) = y0(i-1);

yf(i) = yf(i-1);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%MOVEMENT W/ CB, BEFORE 1st TURN%%%%%%%%%%%%%

if turn_count>0

for i = start/tau:(turn_times(1)/tau-1)

[add, del] = check_for_add(x0(i-1), y0(i-1), xf(i-1), yf(i-1),...

nnx(:, i-1), nny(:, i-1), fil, 0); %check to add

if add == -1

[nnx(:, i), nny(:, i)]= move_wb(turn_angles(1), nnx(:, i-1),...

nny(:, i-1), fil, tau, lam_size, x0(i-1), y0(i-1), xf(i-1), yf(i-1));

edx = sum(nnx(:, i)-nnx(:, i-1))/fil;

x0(i) = x0(i-1)+tau;

y0(i) = 1;

xf(i) = xf(i-1)+tau;

yf(i) = -1;

else

[nnx(:, i), nny(:, i)]= add_fil(nnx(:, i-1), nny(:, i-1),...

x0(i-1), y0(i-1), xf(i-1), yf(i-1), add, del);

x0(i) = x0(i-1);

xf(i) = xf(i-1);

y0(i) = y0(i-1);

yf(i) = yf(i-1);

fprintf('\nI have added a filament number %f after %f timestep\n',add, i);

end

[nnx(:, turn_times(1)/tau), nny(:, turn_times(1)/tau)] = move_wb(turn_angles(1),...

nnx(:, turn_times(1)/tau-1), nny(:, turn_times(1)/tau-1),...

fil, tau, lam_size, x0(turn_times(1)/tau-1),...

y0(turn_times(1)/tau-1), xf(turn_times(1)/tau-1), yf(turn_times(1)/tau-1));

x0(turn_times(1)/tau) = (x0(turn_times(1)/tau-1)+...

xf(turn_times(1)/tau-1))/2-sin(turn_angles(1));

y0(turn_times(1)/tau) = (y0(turn_times(1)/tau-1)+...

yf(turn_times(1)/tau-1))/2+cos(turn_angles(1));

xf(turn_times(1)/tau) = (xf(turn_times(1)/tau-1)+...

x0(turn_times(1)/tau-1))/2+sin(turn_angles(1));

yf(turn_times(1)/tau) = (yf(turn_times(1)/tau-1)+...

y0(turn_times(1)/tau-1))/2-cos(turn_angles(1));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%ALL OTHER MOVEMENT%%%%%%%%%%%%%%%%%%%%%%%

if turn_count == 1

for i = (turn_times(1)/tau+1):time/tau

[add, del] = check_for_add(x0(i-1), y0(i-1), xf(i-1), yf(i-1), nnx(:, i-1), nny(:, i-1),...

fil, turn_angles(1)); %add can equal to 1 through fil-1

%-1-nothing is added, add>0- filament is added between add and add +1

if add == -1

[nnx(:, i), nny(:, i)]= move_wb(turn_angles(1), nnx(:, i-1),...

nny(:, i-1), fil, tau, lam_size, x0(i-1), y0(i-1), xf(i-1), yf(i-1));

edx = sum(nnx(:, i)-nnx(:, i-1))/fil;

edy = sum(nny(:, i)-nny(:, i-1))/fil;

x0(i) = x0(i-1)+cos(turn_angles(1))*tau*.7;

y0(i) = y0(i-1)+sin(turn_angles(1))*tau*.7;

xf(i) = xf(i-1)+cos(turn_angles(1))*tau*.7;

yf(i) = yf(i-1)+sin(turn_angles(1))*tau*.7;

else

[nnx(:, i), nny(:, i)]= add_fil(nnx(:, i-1), nny(:, i-1),...

x0(i-1), y0(i-1), xf(i-1), yf(i-1), add, del);

fprintf('\nI have added a filament number %f after %f timestep\n',add, i);

x0(i) = x0(i-1);

xf(i) = xf(i-1);

y0(i) = y0(i-1);

yf(i) = yf(i-1);

end

elseif turn_count>1

for j = 2:turn_count

for i = turn_times(j-1)/tau:(turn_times(j)/tau-1)

[add, del] = check_for_add(x0(i-1), y0(i-1), xf(i-1), yf(i-1), nnx(:, i-1), nny(:, i-1), fil, sum(turn_angles(1:j)));

if add == -1

[nnx(:, i), nny(:, i)]= move_wb(sum(turn_angles(1:j)),...

nnx(:, i-1), nny(:, i-1), fil, tau, lam_size, x0(i-1), y0(i-1), xf(i-1), yf(i-1));

x0(i) = x0(i-1)+cos(turn_angles(j))*tau*.7;

y0(i) = y0(i-1)+sin(turn_angles(j))*tau*.7;

xf(i) = xf(i-1)+cos(turn_angles(j))*tau*.7;

yf(i) = yf(i-1)+sin(turn_angles(j))*tau*.7;

else

[nnx(:, i), nny(:, i)]= add_fil(nnx(:, i-1), nny(:, i-1), x0(i-1), y0(i-1), xf(i-1), yf(i-1), add, del);

fprintf('\nI have added a filament number %f after %f timestep\n',add, i);

x0(i) = x0(i-1);

xf(i) = xf(i-1);

y0(i) = y0(i-1);

yf(i) = yf(i-1);

end

%%%%%%%%%%%%%%%THE REST OF THE%%%%%%%%%

%%%%%%%%%%%%%%%MOVEMENT%%%%%%%%%%%%%%%%

if (turn_count>0)

for i = turn_times(turn_count)/tau:time/tau

[add, del] = check_for_add(x0(i-1), y0(i-1), xf(i-1), yf(i-1), nnx(:, i-1), nny(:, i-1)...

, fil, sum(turn_angles(1:turn_count)) ); %add can equal to 0 through fil-1

%add=0-nothing is added, add>0- filament is added between add and add +1

if add == -1

[nnx(:, i), nny(:, i)]= move_wb(sum(turn_angles(1:turn_count)),...

nnx(:, i-1), nny(:, i-1), fil, tau, lam_size, x0(i-1), y0(i-1), xf(i-1), yf(i-1));

x0(i) = x0(i-1)+cos(turn_angles(turn_count))*tau*.7;

y0(i) = y0(i-1)+sin(turn_angles(turn_count))*tau*.7;

xf(i) = xf(i-1)+cos(turn_angles(turn_count))*tau*.7;

yf(i) = yf(i-1)+sin(turn_angles(turn_count))*tau*.7;

else

[nnx(:, i), nny(:, i)]= add_fil(nnx(:, i-1), nny(:, i-1), x0(i-1), y0(i-1), xf(i-1), yf(i-1), add, del);

fprintf('\nI have added a filament number %f after %f timestep\n',add, i);

x0(i) = x0(i-1);

xf(i) = xf(i-1);

y0(i) = y0(i-1);

yf(i) = yf(i-1);

end

%m = moviein(time/tau);

%for t =1:time/tau

%plot(nnx(:, 1:t), nny(:, 1:t), '+');

%hold on

%plot([x0(t), xf(t)], [y0(t), yf(t)]);

%hold off

%axis([-.2 10 -3 3]);

% m (:,t) = getframe;

%end

xplot(fil,time/tau) = zeros;

yplot(fil,time/tau) = zeros;

m = moviein(time/tau);

for t =1:time/tau

[xshape,yshape,angle1] = body(x0(t), y0(t), xf(t) , yf(t) , lam_size);

for i = 1: fil

for j = 1:t

if angle1 == 0

if nnx(i,j) <= x0(t)

xplot(i,j) = nnx(i,j);

yplot(i,j) = nny(i,j);

else

xplot(i,j) = 0;

yplot(i,j) = 0;

end

elseif angle1 >0

if nny(i,j) <= ( (y0(t) -yf(t) )/(x0(t) - xf(t) )) *(nnx(i,j) - x0(t) ) +y0(t);

xplot(i,j) = nnx(i,j);

yplot(i,j) = nny(i,j);

else

xplot(i,j) = 0;

yplot(i,j) = 0;

end

else

if nny(i,j) >= ( (y0(t) -yf(t) )/(x0(t) - xf(t) )) *(nnx(i,j) - x0(t) ) +y0(t);

xplot(i,j) = nnx(i,j);

yplot(i,j) = nny(i,j);

else

xplot(i,j) = 0;

yplot(i,j) = 0;

end

plot(nnx(:,1:t),nny(:,1:t),'k.',xplot(:,1:t),yplot(:,1:t),'w.')

axis ([-.2 5 -3 3])

hold on

fill(xshape,yshape,'k');

hold off

m (:,t) = getframe;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subroutine to add filaments… SAVE AS “add_fil.m”

%add_fil.m

%Version 7/23/08

%Brian, Debbie, Kaitlyn

%Subroutine that adds a new filament to the lamellipod and deletes an

%existing filament

%The basic structure of this function is as follows:

%We are provided with an initial vector of x coordinates, which we call xvec_prev, and

%y coordinates which we can yvec_prev, for the initial lamellipod contour. The output

%of the function is a new vector of x coordinates, which we call xvec_next,

%and y coordinates, which we call yvec_next, which are the coordinates for

%the new lamellipod contour, with the added filament and the deleted filament omitted .

%The program code is broken down as follows: We first determine which case

%to use, 1, 2, or 3, then based on which case we need to use, we work through

%steps A-B or A.

%1. We interpolate if the filament removed is NOT THE STARTING OR ENDING

%point of the parabola.

%A. We determine if the filament is on the top half of the parabola

%above the x-axis or the bottom half of the parabola below the x-axis

%B. We chooose linear or parabolic interpolation based on the determinant

%of matrix A and the slopes of the points used for interpolation.

%2. We interpolate if the filament is the STARTING point of the parabola

%A. We chooose linear or parabolic interpolation based on the determinant

%of matrix A and the slopes of the points used for interpolation.

%3. We interpolate if the filament is the ENDING point of the parabola

%A. We chooose linear or parabolic interpolation based on the determinant

%of matrix A and the slopes of the points used for interpolation.

% xvec_prev, yvec_prev are the coordinates of the initial lamellipod contour, x0, y0,

% xf, and yf are the starting and ending points of the lamellipod respectively, and ins is the

% number of the filament that we will take into consideration.

function [xvec_next,yvec_next]= add_fil(xvec_prev, yvec_prev,x0,y0,xf,yf,ins,del)

%n is the number of filaments along the lamellipod

n = length (xvec_prev);

%Case 1: Interpolate the parabola or line if the filament to be removed is not the

%starting or ending point of the parabola.

if ins ~= 0 && ins ~= n

%A. Calculate matrices A and X that will be used for interpolation if

%filament lies on the half of the parabola above the x-axis

if ins < n/2

A = [ yvec_prev(ins)^2, yvec_prev(ins), 1 ; yvec_prev(ins+1)^2, yvec_prev(ins+1), 1 ;...

yvec_prev(ins+2)^2, yvec_prev(ins+2), 1];

X = [xvec_prev(ins) ; xvec_prev(ins+1) ; xvec_prev(ins+2)];

%Define 3 points used for interpolation

y1 = yvec_prev(ins);

y2 = yvec_prev(ins+1);

y3 = yvec_prev(ins+2);

x1 = xvec_prev(ins);

x2 = xvec_prev(ins+1);

x3 = xvec_prev(ins+2);

%Calculate A and X if the filament lies on the half of the parabola below the x-axis

else

A = [ yvec_prev(ins+1)^2, yvec_prev(ins+1), 1 ; yvec_prev(ins)^2, yvec_prev(ins), 1 ;...

yvec_prev(ins-1)^2, yvec_prev(ins-1), 1];

X = [xvec_prev(ins+1) ; xvec_prev(ins) ; xvec_prev(ins-1)];

%Define 3 points used for interpolation

y1 = yvec_prev(ins+1);

y2 = yvec_prev(ins);

y3 = yvec_prev(ins-1);

x1 = xvec_prev(ins+1);

x2 = xvec_prev(ins);

x3 = xvec_prev(ins-1);

end

%B. The determinant of the matrix A that is used to determine the

%coefficents of the parabola must not be (approaching) zero (otherwise it is not

%invertible). Also, the slope of the three points used to

%interpolate the parabola must not be the same.

%If the determinant of A = 0 or the slopes of the three points are

%the same, we use LINEAR interpolation

if abs(det(A))< .0001 || abs ((y2 - y1) / (x2 - x1) - (y3 - y2 ) / (x3 - x2)) < .0001 || isnan(det(A)) == 1

%New x and y points

xstar = (x1 + x2) /2;

ystar = ((xstar - x1)*(y2 - y1))/(x2 - x1) + y1;

if del == 1

xvec_next = [xvec_prev(2:ins) ; xstar ; xvec_prev(ins+1:n)];

yvec_next = [yvec_prev(2:ins) ; ystar ; yvec_prev(ins+1:n)];

else

xvec_next = [xvec_prev(1:ins) ; xstar ; xvec_prev(ins+1:n-1)];

yvec_next = [yvec_prev(1:ins) ; ystar ; yvec_prev(ins+1:n-1)];

end

%If the determinant of A does not equal zero, and the slopes of the

%points are not equal, we use PARABOLIC interpolation

else

B = inv(A);

C = B*X;

%Calculate new x and y points. xstar and ystar are the

%coordinates for the new filament that is added

ystar = (y1+y2) / 2; %New y coordinate

xstar = C(1)*(ystar)^2 + C(2)*ystar + C(3); %New x coordinate

if del == 1

xvec_next = [xvec_prev(2:ins) ; xstar ; xvec_prev(ins+1:n)];

yvec_next = [yvec_prev(2:ins) ; ystar ; yvec_prev(ins+1:n)];

else

xvec_next = [xvec_prev(1:ins) ; xstar ; xvec_prev(ins+1:n-1)];

yvec_next = [yvec_prev(1:ins) ; ystar ; yvec_prev(ins+1:n-1)];

end

%Case 2: Interpolate the parabola or line if the filament to be removed is the

%starting point of the parabola.

elseif ins == 0

if abs (yvec_prev(1) - yvec_prev(2) ) < 0.01

xstar = (x0 + xvec_prev(1) ) /2;

ystar = (y0 + yvec_prev(1) ) /2;

xvec_next = [xvec_prev(1:n-1); xstar];

yvec_next = [yvec_prev(1:n-1); ystar];

else

G = [ yvec_prev(1)*( (xvec_prev(2) - xvec_prev(1))/(yvec_prev(2) - yvec_prev(1))) + xvec_prev(1) ;

y0*( (xf - x0)/(yf - y0)) + x0 ];

H = [ (xvec_prev(2) - xvec_prev(1))/(yvec_prev(2)-yvec_prev(1)) , -1 ; (xf - x0)/(yf - y0) , -1 ];

Z = inv(H)*G;

ystar = (Z(1) + yvec_prev(n))/2;

xstar = (Z(2) + xvec_prev(n))/2;

xvec_next = [xvec_prev(1:n-1); xstar];

yvec_next = [yvec_prev(1:n-1); ystar];

end

%Case 3: Interpolate the parabola or line if the filament to be removed is the

%ending point of the parabola.

else

if abs (yvec_prev(n-1) - yvec_prev(n) ) <0.01

xstar = (xf + xvec_prev(n) ) /2 ;

ystar = (yf + yvec_prev(n) ) /2 ;

xvec_next = [ xvec_prev(2:n) ; xstar];

yvec_next = [ yvec_prev(2:n) ; ystar];

else

G = [ y0 * ( (xf - x0)/(yf - y0) ) + x0 ;

yvec_prev(n) * ( (xvec_prev(n-1) - xvec_prev(n))/(yvec_prev(n-1) - yvec_prev(n)))-xvec_prev(n)];

H = [ (xf - x0) / (yf - y0) , -1 ; (xvec_prev(n-1)-xvec_prev(n)) / (yvec_prev(n-1) - yvec_prev(n)) , -1 ];

Z = inv(H)*G;

ystar = (Z(1) + yvec_prev(n))/2;

xstar = (Z(2) + xvec_prev(n))/2;

xvec_next = [ xvec_prev(2:n) ; xstar];

yvec_next = [ yvec_prev(2:n) ; ystar];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subroutine to describe shape of body-lammelipod interface

%%%%%%%%%%%%%%SAVE AS “body.m”

%body.m

function [xshape,yshape,angle1]= body(x0, y0, xf , yf , lam_size)

yinitial = [ .9 ; .8 ; .7 ; .6 ; .5 ; .4 ; .3 ; .2 ; .1 ; 0 ; -.1 ; -.2 ; -.3 ; -.4 ; -.5 ; -.6 ; -.7 ; - .8 ; -.9 ];

xinitial = (lam_size)*sqrt(1-yinitial.^2);

if xf - x0 == 0

angle1 = 0;

else

angle1 = atan( (y0-yf) / (xf - x0) );

end

distance = xinitial;

if (x0 - xf) == 0

xshape = [x0 ; x0 + distance ; xf ;x0];

yshape = [y0 ; yinitial ; yf; y0];

elseif angle1 >0

m = (y0 - yf) / (x0-xf) ;

ydistance = abs(y0 - yf);

yrear = [y0-ydistance*(1/20) ; y0 - (2/20)*ydistance ; y0 - (3/20)*ydistance ;y0 - (4/20)*ydistance ;y0 - (5/20)*ydistance ;y0 - (6/20)*ydistance ;y0 - (7/20)*ydistance ;y0 - (8/20)*ydistance ;y0 - (9/20)*ydistance ;...

y0 - (10/20)*ydistance ;y0 - (11/20)*ydistance ;y0 - (12/20)*ydistance ;y0 - (13/20)*ydistance ;y0 - (14/20)*ydistance ;y0 - (15/20)*ydistance ;y0 - (16/20)*ydistance ;y0 - (17/20)*ydistance ;y0 - (18/20)*ydistance ;y0 - (19/20)*ydistance ];

xrear = (yrear - y0 +m*x0)/m;

xshape = [x0 ; xrear + (cos(pi/2 - angle1))*distance; xf ; x0] ;

yshape = [y0 ; yrear + sin(pi/2 - angle1)*distance; yf; y0] ;

else

m = (y0 - yf) / (x0-xf) ;

ydistance = abs(y0 - yf);

xrear = (yrear - y0 +m*x0)/m;

xshape = [x0 ; xrear + cos(pi/2 - abs(angle1))*distance; xf; x0 ] ;

yshape = [y0 ; yrear - sin(pi/2 - abs(angle1))*distance; yf; y0 ] ;

end

%%%%%%%%subroutine to check if new bundles should be added based on conservation of momentum

%%%%%%%%%%SAVE AS “check_for_add.m”

% check_for_add.m

function [newfil, delfil]=check_for_add(x0, y0, xf, yf, xvec, yvec, fil, theta)

%Look at first point

%Get equation of the line

%Check is points are to one or opposite sides

%if opposite check for greatest distance

%Look at last

%blah blah

%fprintf('\ny0 = %f, x0 = %f\n', length(y0), length(x0));

delfil=-1;

newfil=-1;

while(1 == 1)

%EQUATION OF THE LINE

A = (yf-y0)/(xf-x0); %A is the slope

B = (y0-A*x0); % intersection w/ y axis

if xf == x0

if (xvec(1)-x0)*(xvec(2)-x0)<0

if xvec(1)==xvec(fil)

dist_vect =((xvec(3:fil-1)-xvec(2:fil-2)).^2+(yvec(3:fil-1)-yvec(2:fil-2)).^2).^(1/2);

[val, coord] = max(dist_vect);

newfil = coord+1;

delfil = 1;

break

else

dist_vect =((xvec(3:fil)-xvec(2:fil-1)).^2+(yvec(3:fil)-yvec(2:fil-1)).^2).^(1/2);

[val, coord] = max(dist_vect);

newfil = coord+1;

delfil = 1;

break

end

if (xvec(fil)-x0)*(xvec(fil-1)-x0)<0

dist_vect =((xvec(1:fil-2)-xvec(2:fil-1)).^2+(yvec(1:fil-2)-yvec(2:fil-1)).^2).^(1/2);

[val, coord] = max(dist_vect);

newfil = coord;

delfil = fil;

break

end

else

if (yvec(1)-A*xvec(1)-B)*(yvec(2)-A*xvec(2)-B)<0

dist_vect =((xvec(3:fil)-xvec(2:fil-1)).^2+(yvec(3:fil)-yvec(2:fil-1)).^2).^(1/2);

[val, coord] = max(dist_vect);

newfil = coord+1;

delfil = 1;

break

end

if (yvec(fil)-A*xvec(fil)-B)*(yvec(fil-1)-A*xvec(fil-1)-B)<0

dist_vect =((xvec(1:fil-2)-xvec(2:fil-1)).^2+(yvec(1:fil-2)-yvec(2:fil-1)).^2).^(1/2);

[val, coord] = max(dist_vect);

newfil = coord;

delfil = fil;

break

end

newfil = -1;

delfil = -1;

break

end

%%%%%%%%%%%%%%%%%%%%%%%subroutine to move the lamellipod and body SAVE as “move_wb.m”

% move_wb.m

function[xvec_next, yvec_next]= move_wb(angle, xvec_prev, yvec_prev, fil, tau, lam_size, x0, y0, xf, yf)

%%%%%THIS FUNCTION LOOKS AT THE MOTION OF THE LAMELLIPOD AFTER THE CELL

%%%%%BODY MOTION IS TRIGGERED%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if length(fil)~=1 || length(angle)~=1 || length(tau)~=1 || length(lam_size)~=1

error('Error #1'); %what needs to be a scalar is a vector

end

if length(xvec_prev)~=fil || length(yvec_prev)~=fil

error('Error #2'); %vector dimentions are incorrect

end

if (angle>pi/2 || angle<-pi/2)

error('Error #3'); %too much turning

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xvec_previous = [xvec_prev(1); xvec_prev; xvec_prev(fil)];

yvec_previous = [yvec_prev(1); yvec_prev; yvec_prev(fil)];

if (prod(yvec_previous(1:fil)- yvec_previous(3:(fil+2)))== 0)

error('Error # 4'); %dividing by zero somewhere

end

slope = -(xvec_previous(1:(fil))- xvec_previous(3:(fil+2)))./...

(yvec_previous(1:fil)- yvec_previous(3:(fil+2)));

cosn = 1./(sqrt(1+slope.^2));

sinn = slope./(sqrt(1+slope.^2));

pH_vec = pH_function2(angle, xvec_prev, yvec_prev, fil, x0, y0, xf, yf);

xvec_next = xvec_prev+tau*cosn.* pH_vec;

yvec_next = yvec_prev+tau*sinn.* pH_vec;

for i = 1:(fil-1)/2

if yvec_next(i)<yvec_prev(i)

yvec_next(i) = yvec_prev(i);

xvec_next(i) = xvec_prev(i);

end

for i = (fil+1)/2:fil

if yvec_next(i)>yvec_prev(i)

yvec_next(i) = yvec_prev(i);

xvec_next(i) = xvec_prev(i);

end

%%%%%%%%%%%%%%%%%%%%%%%%subroutine to extend the lamellipod SAVE AS “move_wob.m”

%move_wob.m

function[xvec_next, yvec_next]= move_wob(xvec_prev, yvec_prev, fil, tau, lam_size)

%we only have the growing filaments be part of the vector

% THE WAY THIS FUNCTION IS DIFFERENT FROM WOB IS THAT IT ASSUMES THE

% EXTENTION AND ELONGATION OF THE LAMELLIPOD THAT EQUALS TO THE SPEED OF

% MOTION OF THE FRONT FILAMENT%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%CHECK IF ALL THE DIMENTIONS ARE CORRECT%%%%%

if (length(fil)~=1 ||length(tau)~=1 ||length(lam_size)~=1)

error('Error #1'); %what needs to be a scalar is a vector

end

if (length(xvec_prev)~=fil || length(yvec_prev)~=fil)

error('Error #2'); %vector dimentions are incorrect

end

%%%%AT THIS MOMENT WE ASSUME THAT THE TURNING HAPPENS ONLY AFTER THE CELL

%%%%BODY IS TRIGGERED TO MOVE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%adding two points to the vectors at the front and back to calculate slopes

%at all the points

%%%MAYBE IT IS BETTER TO USE A BACKWARD AND A FORWARD DIFFERENCE AT THOSE

%%%POINTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xvec_previous = [0; xvec_prev; 0];

yvec_previous = [1; yvec_prev; -1];

if (prod(yvec_previous(1:fil)- yvec_previous(3:(fil+2)))== 0)

error('Error # 4'); %dividing by zero somewhere

end

slope = -(xvec_previous(1:(fil))- xvec_previous(3:(fil+2)))./...

(yvec_previous(1:fil)- yvec_previous(3:(fil+2)));

cosn = 1./(sqrt(1+slope.^2)); %cos of the angle of the direction of motion

sinn = slope./(sqrt(1+slope.^2)); %sin of the angle of the direction of motion

for i=1:fil

%fprintf('\ncosn(%f)= %f\n', i, cosn(i));

%fprintf('\nsinn(%f)= %f\n', i, sinn(i));

end

%%%xvec_prev = xvec_prev(2:(fil+1));

%%%yvec_prev = yvec_prev(2:(fil+1));

%pH_vec = pH_function(xvec_prev, yvec_prev, 0, lam_size, fil, x0, y0, xf, yf);

xvec_next = xvec_prev + tau*cosn;%.%* pH_vec; %pH is a separate function

yvec_next = yvec_prev + tau*sinn;%.* pH_vec;

for i = 1:(fil-1)/2

if yvec_next(i)<yvec_prev(i)

yvec_next(i) = yvec_prev(i);

xvec_next(i) = xvec_prev(i);

end

for i = (fil+1)/2:fil

if yvec_next(i)>yvec_prev(i)

yvec_next(i) = yvec_prev(i);

xvec_next(i) = xvec_prev(i);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subroutine for one variation of pH regulation

%%%%%%%%%%%%%%%%%SAVE AS “pH_function.m”

%pH_function.m

%version 7/22/08

%Jerry, Nicole

%Function that determines the magnitude vector based on pH

function [mag_vec]=pH_function(xvec, yvec, theta, lam_size, fil, x0, y0, xf, yf)

%theta is rotation of the front of the cell body

%(x0,y0) is the location of the middle point of the cell body front

%fil is the number of filaments

%Initialize

mag_vec(1:fil) = zeros;

%calculate the middle point of the cell body interface

%x0=xvec((fil+1)/2)-cos(theta)*lam_size;

%y0=yvec((fil+1)/2)-sin(theta)*lam_size;

for i=1:fil

%calculate the distance from cell body interface

%to tip of growing filament

dist=sqrt((xvec(i)-(x0+xf)/2)^2+(yvec(i)-(y0+yf)/2)^2);

%Calculate pH based on distance

pH=6.2*(1+0.03226*dist);

%Determine if pH is in the growth range

%and puts growth into magnitude vector

if pH>=6.30 && pH<=6.4

mag_vec(i)=10*(pH-6.3);

else

mag_vec(i)=0;

end

%%%%%%%%%%%%%%%%%%%%%subroutine for another variation of pH regulation SAVE AS “pH_function1.m”

%pH_function1.m

function [ans]=pH_function1(angle, xvec, yvec, fil, x0, y0, xf, yf)

%Initialize

ans(1:fil) = zeros;

%slope of the body line(A) and the intersection with the y-axis(B)

if(abs(xf-x0)>.00001)

A = (y0-yf)/(x0-xf);

B = y0-A*x0;

%slope of the perpendicular(a1) and the intersection with the y-axis(b1)

a1(1:fil) = -1/A; %row vector

b1(1:fil) = yvec-a1'.*xvec; %column vector

xint = (B-b1)./(a1-A);

xdist = abs(xint'-xvec);

a = xdist./cos(angle);

K = .03;

l = 3;

if a>1/4

ans = 31*(K*a.*exp(-l*(a-1).^2));

else

ans = 0;

end

else

K = .03;

l = 3;

a = abs(x0-xvec);

if a>1/4

ans = 31*(K*a.*exp(-l*(a-1).^2));

else

ans = 0;

end

%%%%%%%%%%%%%%%%%%subroutine for another variation of pH regulation SAVE AS “pH_function2.m”

%pH_function2.m

function [ans]=pH_function2(angle, xvec, yvec, fil, x0, y0, xf, yf, lam_size)

%Initialize

ans(1:fil) = zeros;

%slope of the body line(A) and the intersection with the y-axis(B)

if(abs(xf-x0)>.00001)

A = (y0-yf)/(x0-xf);

B = y0-A*x0;

%slope of the perpendicular(a1) and the intersection with the y-axis(b1)

a1(1:fil) = -1/A; %row vector

b1(1:fil) = yvec-a1'.*xvec; %column vector

xint = (B-b1)./(a1-A);

xdist = abs(xint'-xvec);

ans = xdist/cos(angle);

%ans=3* exp((-.4)*ans).*ans;

ans(floor(fil*.65):fil) = ans(floor(fil*.65):fil)*.1;

else

a = abs(x0-xvec);

ans = 3*exp(-.4*a.^2).*a;

for i = 1:fil

if abs(yvec(i))>1.4

ans(i) = ans(i)*.1;

end

%%%%%%%%%%%%%%%%%%% subroutine to establish initial shape of cell body SAVE AS “shape_init.m”

%shape_init.m

function[xvec_first, yvec_first]= shape_init(choice, coeff, fil)

%%%%%%%%%%%CHOOSING THE INITIAL POINTS USING ANGLE DIVISION%%%%%%%%

%%%(INTERSECTION OF THE ELIPSE WITH LINES GOING THROUGH ZERO)%%%%%%

if choice == 1

for i = 1:fil

slope(i) = tan(((fil+1)/2-i)*pi/(fil+1));

xvec_first(i) = 1/(sqrt((slope(i))^2+coeff^(-2)));

yvec_first(i) = slope(i)*xvec_first(i);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%CHOOSING THE INITIAL POINTS USING Y-AXES DIVISION%%%%%%%

%%%%%(INTERSECTION OF THE ELIPSE WITH LINES PARALLEL TO X-AXIS)%%%

else

for i = 1:fil

yvec_first(i) = -1+ 2*(i/(fil+1));

xvec_first(i) = sqrt((1-yvec_first(i)^2)*(coeff^2));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%subroutine to initiate motion SAVE AS “start_move.m”

% start_move.m

function[t] = start_move(lam_size, tau)

t = .4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%